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

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

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


19
20
21
22
23
24
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
25
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
        plane: 
30
31
    Returns:
        points, direction
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
32
            points (list of arrays of length 3): 3d points
33
34
35
    """
    distance0 = _dist_from_plane(p0, plane)
    distance1 = _dist_from_plane(p1, plane)
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
36
37
    p0_on_plane = abs(distance0) < np.finfo(float).eps
    p1_on_plane = abs(distance1) < np.finfo(float).eps
38
    points = []
39
    direction = 0
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
40
    if p0_on_plane:
41
        points.append(p0)
42

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

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

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

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


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

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

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


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

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

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
141

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

145

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


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

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

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

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

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

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

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

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
        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
314
315
316
317
318
319
320
321
322
323
324
    @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
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
    @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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    def planes(self):
        return self._planes
549

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

553
    def in_faces(self, points, delta, assign_multiple=False):
554
555
        """
        Check whether points lie within triangles with Barycentric Technique
556
        see Triangles3D.in_triangles
557
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
558
        masks = self.triangles().in_triangles(points, delta,
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
559
                                              assign_multiple=assign_multiple)
560
561
        return masks

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
562
    def removeFaces(self, face_delete_mask):
563
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
564
        Remove faces where face_delete_mask is True
565
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
566
567
568
        face_delete_mask = np.array(face_delete_mask, dtype=bool)
        self.faces = self.faces[~face_delete_mask]
        self.faceScalars = self.faceScalars[~face_delete_mask]
569

570
    def template(self, sub_mesh):
571
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
572
        'Manual' way to build a template that can be used with self.cut
573
        Returns:
574
575
            Mesh3D: template (see cut), can be used as template to retrieve
                sub_mesh from self instance
576
        Examples:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
577
578
            >>> mp = tfields.TensorFields([[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
            ...                           [1, 2, 3, 4])
579
            >>> 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
580
            ...                     maps=[mp])
581
            >>> from sympy.abc import y
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
582
583
584
            >>> m_cut = m.cut(y < 1.5, at_intersection='split')
            >>> template = m.template(m_cut)
            >>> assert m_cut.equal(m.cut(template))
585

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
586
587
        TODO:
            fields template not yet implemented
588
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
589
590
        face_indices = np.arange(self.maps[0].shape[0])
        cents = tfields.Tensors(sub_mesh.centroids())
591
        mask = self.in_faces(cents, delta=None)
592
        inst = sub_mesh.copy()
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
593
594
595
596
597
598
599
600
601
602
603
        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)
                        ]
604
605
        return inst

606
    def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False):
607
        """
608
        Partition the mesh with the cuts given and return the template
609
        """
610
611
612
613
614
        eps = 0.000000001
        # direct return if self is empty
        if len(self) == 0:
            return self.copy(), self.copy()

615
616
617
618
619
620
621
622
623
624
625
        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))))

626
        # mask for points that do not fulfill the cut expression
627
        mask = inst.evalf(expression)
628
        # remove the points
629

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

686
                for exprPart in expression_parts:
687
688
689
                    inst, _ = inst._cut_sympy(exprPart,
                                              at_intersection='split',
                                              _in_recursion=True)
690
            elif len(expression_parts) == 1:
691
                # TODO maps[0] -> smthng like inst.get_map(dim=3)
692
693
694
                points = [sympy.symbols('x0, y0, z0'),
                          sympy.symbols('x1, y1, z1'),
                          sympy.symbols('x2, y2, z2')]
695
                plane_sympy = tfields.lib.symbolics.to_plane(expression)
696
697
698
                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}
699

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
700
                norm_vectors = inst.triangles().norms()
701
702
                new_points = np.empty((0, 3))
                new_faces = np.empty((0, 3))
703
                new_fields = [tfields.Tensors(np.empty((0,) + field.shape[1:]),
704
                                              coord_sys=field.coord_sys)
705
706
                              for field in inst.fields]
                new_map_fields = [[] for field in inst.maps[0].fields]
707
                new_norm_vectors = []
708
                newScalarMap = []
709
                n_new = 0
710

711
712
713
714
715
716
717
                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]):
718
                    """
719
                    vertices_rejected is a mask for each face that is True, where
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
720
                    a point is on the rejected side of the plane
721
                    """
722
                    vertices_rejected = [~mask[f] for f in face]
723
724
725
726
727
                    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
728
729
                        n_true = vertices_rejected.count(True)
                        lonely_bool = True if n_true == 1 else False
730

Daniel Boeckenhoff's avatar
dunno    
Daniel Boeckenhoff committed
731
                        triangle_points = [vertices[f] for f in face]
732
                        """
733
                        Add the intersection points and faces
734
                        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
735
                        intersection = _intersect(triangle_points, plane, vertices_rejected)
736
                        last_idx = len(vertices) - 1
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
737
                        for tri_list in intersection:
738
                            new_face = []
739
740
741
                            for item in tri_list:
                                if isinstance(item, int):
                                    # reference to old vertex
742
                                    new_face.append(face[item])
743
744
745
                                elif isinstance(item, complex):
                                    # reference to new vertex that has been
                                    # concatenated already
746
                                    new_face.append(last_idx + int(item.imag))
747
748
                                else:
                                    # new vertex
749
                                    new_face.append(len(vertices))
750
                                    vertices = np.append(vertices,
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
751
                                                         [[float(x) for x in item]],
752
753
754
755
756
                                                         axis=0)
                                    fields = [np.append(field,
                                                        np.full((1,) + field.shape[1:], np.nan),
                                                        axis=0)
                                              for field in fields]
757
                            faces = np.append(faces, [new_face], axis=0)
758
                            faces_fields = [np.append(field,
759
                                                      [field[i]],
760
761
762
763
764
                                                      axis=0)
                                            for field in faces_fields]
                            faces_fields[-1][-1] = i

                face_map = tfields.TensorFields(faces, *faces_fields,
765
                                                dtype=int,
766
                                                coord_sys=inst.maps[0].coord_sys)
767
768
769
                inst = tfields.Mesh3D(vertices,
                                      *fields,
                                      maps=[face_map] + inst.maps[1:],
770
                                      coord_sys=inst.coord_sys)
771
772
773
774
775
                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]
776
            else:
777
                raise ValueError("Sympy expression is not splitable.")
778
            inst = inst.cleaned()
779
        elif at_intersection == 'remove':
780
            inst = inst[mask]
781
        else:
782
783
            raise AttributeError("No at_intersection method called {at_intersection} "
                                 "implemented".format(**locals()))
784
785
786
787
788
789
790
791
792
793
794
795
796

        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)
797
        return inst, template
798
799

    def _cut_template(self, template):
800
801
802
803
804
805
806
807
        """
        Args:
            template (tfields.Mesh3D)

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

808
            Build mesh
809
810
811
            >>> 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],
812
813
            ...                    [0.0, 0.1, 0.2, 0.3, 0.4],
            ...                    [0.0, -0.1, -0.2, -0.3, -0.4],
814
815
            ...                    maps=[mmap])

816
            Build template
817
818
819
            >>> tmap = tfields.TensorFields([[0, 3, 4], [0, 1, 2]],
            ...                             [1, 0])
            >>> t = tfields.Mesh3D([[0]*3, [-1]*3, [-2]*3, [-3]*3, [-4]*3],
820
            ...                    [1, 0, 3, 2, 4],
821
822
            ...                    maps=[tmap])

823
            Use template as instruction to make a fast cut
824
825
826
827
828
829
830
831
832
            >>> 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]])
                                   
        """
833
834
        # Possible Extension (small todo): check: len(field(s)) == len(self/maps)

835
        # Redirect fields
836
        fields = []
837
        if template.fields:
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
            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)
854

855
        # Redirect maps and their fields
856
857
        maps = []
        for mp, template_mp in zip(self.maps, template.maps):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
858
            mp_fields = []
859
860
861
862
            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
863
                    mp_fields.append(field[template_mp.fields[0].astype(int)])
864
            new_mp = tfields.TensorFields(tfields.Tensors(template_mp),
865
866
867
                                          *mp_fields)
            maps.append(new_mp)

868
869
        inst = tfields.Mesh3D(tfields.Tensors(template),
                              *fields,
870
                              maps=maps)
871
872
        return inst

873
    def cut(self, expression, coord_sys=None, at_intersection=None,
874
            return_template=False):
875
876
877
        """
        cut method for Mesh3D.
        Args:
878
879
880
881
882
883
884
885
886
887
            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.
888
            coord_sys (coordinate system to cut in):
889
            at_intersection (str): instruction on what to do, when a cut will intersect a triangle.
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
890
891
892
                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.
893
894
            return_template (bool): If True: return the template
                            to redo the same cut fast
895
896
        Examples:
            define the cut
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
897
            >>> import numpy as np
898
            >>> import tfields
899
            >>> from sympy.abc import x,y,z
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
900
            >>> cut_expr = x > 1.5
901

902
903
            >>> m = tfields.Mesh3D.grid((0, 3, 4),
            ...                         (0, 3, 4),
904
            ...                         (0, 0, 1))
905
906
907
908
909
910
            >>> 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
911
            >>> mNew = m.cut(cut_expr)
912
            >>> len(mNew)
913
            8
914
            >>> mNew.nfaces()
915
916
917
918
            6
            >>> float(mNew[:, 0].min())
            2.0

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
919
920
921
922
923
924
925
926
927
928
            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')
929
            >>> float(m_split[:, 0].min())
930
            1.5
931
            >>> len(m_split)
932
            15
933
            >>> m_split.nfaces()
934
935
            15

936
937
            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
938
            >>> m_split_2, template = m.cut(cut_expr, at_intersection='split',
939
            ...                                    return_template=True)
940
941
942
943
944
945
946
947
948
949
            >>> 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)
950

951
            This seems irrelevant at first but consider, the map field or the
952
953
954
955
956
957
958
959
            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:
960
961
            >>> 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')
962
963
964

        Returns:
            copy of cut mesh
965
            * optional: template
966
967

        """
968
        with self.tmp_transform(coord_sys or self.coord_sys):
969
            if isinstance(expression, Mesh3D):
970
971
                template = expression
                obj = self._cut_template(template)
972
973
974
975
976
977
            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
978

979
980
981
982
983
984
985
986
987
988
989
990
991
    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
992
    def plot(self, **kwargs):  # pragma: no cover
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
993
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
994
        Forwarding to plotTools.plot_mesh
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
995
996
997
998
999
1000
        """
        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]
For faster browsing, not all history is shown. View entire blame