vedo.dolfin

  1#!/usr/bin/env python3
  2# -*- coding: utf-8 -*-
  3import numpy as np
  4
  5try:
  6    import vedo.vtkclasses as vtk
  7except ImportError:
  8    import vtkmodules.all as vtk
  9
 10import vedo
 11from vedo.colors import printc
 12from vedo import utils
 13from vedo import shapes
 14from vedo.mesh import Mesh
 15from vedo.plotter import show
 16
 17__docformat__ = "google"
 18
 19__doc__ = """
 20Submodule for support of the [FEniCS/Dolfin](https://fenicsproject.org) library.
 21
 22Example:
 23    ```python
 24    import dolfin
 25    from vedo import dataurl, download
 26    from vedo.dolfin import plot
 27    fname = download(dataurl+"dolfin_fine.xml")
 28    mesh = dolfin.Mesh(fname)
 29    plot(mesh)
 30    ```
 31
 32![](https://user-images.githubusercontent.com/32848391/53026243-d2d31900-3462-11e9-9dde-518218c241b6.jpg)
 33
 34.. note::
 35    Find many more examples in directory    
 36    [vedo/examples/dolfin](https://github.com/marcomusy/vedo/blob/master/examples/other/dolfin).
 37"""
 38
 39__all__ = [
 40    "plot",
 41]
 42
 43
 44##########################################################################
 45def _inputsort(obj):
 46
 47    import dolfin
 48
 49    u = None
 50    mesh = None
 51    if not utils.is_sequence(obj):
 52        obj = [obj]
 53
 54    for ob in obj:
 55        inputtype = str(type(ob))
 56
 57        # printc('inputtype is', inputtype, c=2)
 58
 59        if "vedo" in inputtype:  # skip vtk objects, will be added later
 60            continue
 61
 62        if "dolfin" in inputtype or "ufl" in inputtype:
 63            if "MeshFunction" in inputtype:
 64                mesh = ob.mesh()
 65
 66                if ob.dim()>0:
 67                    printc('MeshFunction of dim>0 not supported.', c='r')
 68                    printc('Try e.g.:  MeshFunction("size_t", mesh, 0)', c='r', italic=1)
 69                    printc('instead of MeshFunction("size_t", mesh, 1)', c='r', strike=1)
 70                else:
 71                    # printc(ob.dim(), mesh.num_cells(), len(mesh.coordinates()), len(ob.array()))
 72                    V = dolfin.FunctionSpace(mesh, "CG", 1)
 73                    u = dolfin.Function(V)
 74                    v2d = dolfin.vertex_to_dof_map(V)
 75                    u.vector()[v2d] = ob.array()
 76            elif "Function" in inputtype or "Expression" in inputtype:
 77                u = ob
 78            elif "ufl.mathfunctions" in inputtype:  # not working
 79                u = ob
 80            elif "Mesh" in inputtype:
 81                mesh = ob
 82            elif "algebra" in inputtype:
 83                mesh = ob.ufl_domain()
 84                # print('algebra', ob.ufl_domain())
 85
 86    if u and not mesh and hasattr(u, "function_space"):
 87        V = u.function_space()
 88        if V:
 89            mesh = V.mesh()
 90    if u and not mesh and hasattr(u, "mesh"):
 91        mesh = u.mesh()
 92
 93    # printc('------------------------------------')
 94    # printc('mesh.topology dim=', mesh.topology().dim())
 95    # printc('mesh.geometry dim=', mesh.geometry().dim())
 96    # if u: printc('u.value_rank()', u.value_rank())
 97    # if u and u.value_rank(): printc('u.value_dimension()', u.value_dimension(0)) # axis=0
 98    ##if u: printc('u.value_shape()', u.value_shape())
 99    return (mesh, u)
100
101
102def _compute_uvalues(u, mesh):
103    # the whole purpose of this function is
104    # to have a scalar (or vector) for each point of the mesh
105
106    if not u:
107        return np.array([])
108    #    print('u',u)
109
110    if hasattr(u, "compute_vertex_values"):  # old dolfin, works fine
111        u_values = u.compute_vertex_values(mesh)
112
113        if u.value_rank() and u.value_dimension(0) > 1:
114            l = u_values.shape[0]
115            u_values = u_values.reshape(u.value_dimension(0), int(l/u.value_dimension(0))).T
116
117    elif hasattr(u, "compute_point_values"):  # dolfinx
118        u_values = u.compute_point_values()
119
120        try:
121            from dolfin import fem
122
123            fvec = u.vector
124        except RuntimeError:
125            fspace = u.function_space
126            try:
127                fspace = fspace.collapse()
128            except RuntimeError:
129                return []
130
131        fvec = fem.interpolate(u, fspace).vector
132
133        tdim = mesh.topology.dim
134
135        # print('fvec.getSize', fvec.getSize(), mesh.num_entities(tdim))
136        if fvec.getSize() == mesh.num_entities(tdim):
137            # DG0 cellwise function
138            C = fvec.get_local()
139            if C.dtype.type is np.complex128:
140                print("Plotting real part of complex data")
141                C = np.real(C)
142
143        u_values = C
144
145    else:
146        u_values = []
147
148    if hasattr(mesh, "coordinates"):
149        coords = mesh.coordinates()
150    else:
151        coords = mesh.geometry.points
152
153    if u_values.shape[0] != coords.shape[0]:
154        vedo.logger.warning("mismatch in vedo.dolfin._compute_uvalues")
155        u_values = np.array([u(p) for p in coords])
156    return u_values
157
158
159def plot(*inputobj, **options):
160    """
161    Plot the object(s) provided.
162
163    Input can be any combination of: `Mesh`, `Volume`, `dolfin.Mesh`,
164    `dolfin.MeshFunction`, `dolfin.Expression` or `dolfin.Function`.
165
166    Return the current `Plotter` class instance.
167
168    Arguments:
169        mode : (str)
170            one or more of the following can be combined in any order
171            - `mesh`/`color`, will plot the mesh, by default colored with a scalar if available
172            - `displacement` show displaced mesh by solution
173            - `arrows`, mesh displacements are plotted as scaled arrows.
174            - `lines`, mesh displacements are plotted as scaled lines.
175            - `tensors`, to be implemented
176        add : (bool)
177            add the input objects without clearing the already plotted ones
178        density : (float)
179            show only a subset of lines or arrows [0-1]
180        wire[frame] : (bool)
181            visualize mesh as wireframe [False]
182        c[olor] : (color)
183            set mesh color [None]
184        exterior : (bool)
185            only show the outer surface of the mesh [False]
186        alpha : (float)
187            set object's transparency [1]
188        lw : (int)
189            line width of the mesh (set to zero to hide mesh) [0.1]
190        ps :  int
191            set point size of mesh vertices [None]
192        z : (float)
193            add a constant to z-coordinate (useful to show 2D slices as function of time)
194        legend : (str)
195            add a legend to the top-right of window [None]
196        scalarbar : (bool)
197            add a scalarbar to the window ['vertical']
198        vmin : (float)
199            set the minimum for the range of the scalar [None]
200        vmax : (float)
201            set the maximum for the range of the scalar [None]
202        scale : (float)
203            add a scaling factor to arrows and lines sizes [1]
204        cmap : (str)
205            choose a color map for scalars
206        shading : (str)
207            mesh shading ['flat', 'phong']
208        text : (str)
209            add a gray text comment to the top-left of the window [None]
210        isolines : (dict)
211            dictionary of isolines properties
212            - n, (int) - add this number of isolines to the mesh
213            - c, - isoline color
214            - lw, (float) - isoline width
215            - z, (float) - add to the isoline z coordinate to make them more visible
216        streamlines : (dict)
217            dictionary of streamlines properties
218            - probes, (list, None) - custom list of points to use as seeds
219            - tol, (float) - tolerance to reduce the number of seed points used in mesh
220            - lw, (float) - line width of the streamline
221            - direction, (str) - direction of integration ('forward', 'backward' or 'both')
222            - max_propagation, (float) - max propagation of the streamline
223            - scalar_range, (list) - scalar range of coloring
224        warpZfactor : (float)
225            elevate z-axis by scalar value (useful for 2D geometries)
226        warpYfactor : (float)
227            elevate z-axis by scalar value (useful for 1D geometries)
228        scaleMeshFactors : (list)
229            rescale mesh by these factors [1,1,1]
230        new : (bool)
231            spawn a new instance of Plotter class, pops up a new window
232        at : (int)
233            renderer number to plot to
234        shape : (list)
235            subdvide window in (n,m) rows and columns
236        N : (int)
237            automatically subdvide window in N renderers
238        pos : (list)
239            (x,y) coordinates of the window position on screen
240        size : (list)
241            window size (x,y)
242        title : (str)
243            window title
244        bg : (color)
245            background color name of window
246        bg2 : (color)
247            second background color name to create a color gradient
248        style : (int)
249            choose a predefined style [0-4]
250            - 0, `vedo`, style (blackboard background, rainbow color map)
251            - 1, `matplotlib`, style (white background, viridis color map)
252            - 2, `paraview`, style
253            - 3, `meshlab`, style
254            - 4, `bw`, black and white style.
255        axes : (int)
256            Axes type number.
257            Axes type-1 can be fully customized by passing a dictionary `axes=dict()`.
258            - 0,  no axes,
259            - 1,  draw customizable grid axes (see below).
260            - 2,  show cartesian axes from (0,0,0)
261            - 3,  show positive range of cartesian axes from (0,0,0)
262            - 4,  show a triad at bottom left
263            - 5,  show a cube at bottom left
264            - 6,  mark the corners of the bounding box
265            - 7,  draw a simple ruler at the bottom of the window
266            - 8,  show the `vtkCubeAxesActor` object,
267            - 9,  show the bounding box outLine,
268            - 10, show three circles representing the maximum bounding box,
269            - 11, show a large grid on the x-y plane (use with zoom=8)
270            - 12, show polar axes.
271        infinity : (bool)
272            if True fugue point is set at infinity (no perspective effects)
273        sharecam : (bool)
274            if False each renderer will have an independent vtkCamera
275        interactive : (bool)
276            if True will stop after show() to allow interaction w/ window
277        offscreen : (bool)
278            if True will not show the rendering window
279        zoom : (float)
280            camera zooming factor
281        viewup : (list), str
282            camera view-up direction ['x','y','z', or a vector direction]
283        azimuth : (float)
284            add azimuth rotation of the scene, in degrees
285        elevation : (float)
286            add elevation rotation of the scene, in degrees
287        roll : (float)
288            add roll-type rotation of the scene, in degrees
289        camera : (dict)
290            Camera parameters can further be specified with a dictionary
291            assigned to the `camera` keyword:
292            (E.g. `show(camera={'pos':(1,2,3), 'thickness':1000,})`)
293            - `pos`, `(list)`,
294                the position of the camera in world coordinates
295            - `focal_point`, `(list)`,
296                the focal point of the camera in world coordinates
297            - `viewup`, `(list)`,
298                the view up direction for the camera
299            - `distance`, `(float)`,
300                set the focal point to the specified distance from the camera position.
301            - `clipping_range`, `(float)`,
302                distance of the near and far clipping planes along the direction of projection.
303            - `parallel_scale`, `(float)`,
304                scaling used for a parallel projection, i.e. the height of the viewport
305                in world-coordinate distances. The default is 1. Note that the "scale" parameter works as
306                an "inverse scale", larger numbers produce smaller images.
307                This method has no effect in perspective projection mode.
308            - `thickness`, `(float)`,
309                set the distance between clipping planes. This method adjusts the far clipping
310                plane to be set a distance 'thickness' beyond the near clipping plane.
311            - `view_angle`, `(float)`,
312                the camera view angle, which is the angular height of the camera view
313                measured in degrees. The default angle is 30 degrees.
314                This method has no effect in parallel projection mode.
315                The formula for setting the angle up for perfect perspective viewing is:
316                angle = 2*atan((h/2)/d) where h is the height of the RenderWindow
317                (measured by holding a ruler up to your screen) and d is the distance
318                from your eyes to the screen.
319    """
320    if len(inputobj) == 0:
321        vedo.plotter_instance.interactive()
322        return None
323
324    if "numpy" in str(type(inputobj[0])):
325        from vedo.pyplot import plot as pyplot_plot
326
327        return pyplot_plot(*inputobj, **options)
328
329    mesh, u = _inputsort(inputobj)
330
331    mode = options.pop("mode", "mesh")
332    ttime = options.pop("z", None)
333
334    add = options.pop("add", False)
335
336    wire = options.pop("wireframe", None)
337
338    c = options.pop("c", None)
339    color = options.pop("color", None)
340    if color is not None:
341        c = color
342
343    lc = options.pop("lc", None)
344
345    alpha = options.pop("alpha", 1)
346    lw = options.pop("lw", 0.5)
347    ps = options.pop("ps", None)
348    legend = options.pop("legend", None)
349    scbar = options.pop("scalarbar", "v")
350    vmin = options.pop("vmin", None)
351    vmax = options.pop("vmax", None)
352    cmap = options.pop("cmap", None)
353    scale = options.pop("scale", 1)
354    scaleMeshFactors = options.pop("scaleMeshFactors", [1, 1, 1])
355    shading = options.pop("shading", "phong")
356    text = options.pop("text", None)
357    style = options.pop("style", "vtk")
358    isolns = options.pop("isolines", {})
359    streamlines = options.pop("streamlines", {})
360    warpZfactor = options.pop("warpZfactor", None)
361    warpYfactor = options.pop("warpYfactor", None)
362    lighting = options.pop("lighting", None)
363    exterior = options.pop("exterior", False)
364    returnActorsNoShow = options.pop("returnActorsNoShow", False)
365    at = options.pop("at", 0)
366
367    # refresh axes titles for axes type = 8 (vtkCubeAxesActor)
368    xtitle = options.pop("xtitle", "x")
369    ytitle = options.pop("ytitle", "y")
370    ztitle = options.pop("ztitle", "z")
371    if vedo.plotter_instance:
372        if xtitle != "x":
373            aet = vedo.plotter_instance.axes_instances
374            if len(aet) > at and isinstance(aet[at], vtk.vtkCubeAxesActor):
375                aet[at].SetXTitle(xtitle)
376        if ytitle != "y":
377            aet = vedo.plotter_instance.axes_instances
378            if len(aet) > at and isinstance(aet[at], vtk.vtkCubeAxesActor):
379                aet[at].SetYTitle(ytitle)
380        if ztitle != "z":
381            aet = vedo.plotter_instance.axes_instances
382            if len(aet) > at and isinstance(aet[at], vtk.vtkCubeAxesActor):
383                aet[at].SetZTitle(ztitle)
384
385    # change some default to emulate standard behaviours
386    if style in (0, 'vtk'):
387        axes = options.pop("axes", None)
388        if axes is None:
389            options["axes"] = {
390                "xygrid": False,
391                "yzgrid": False,
392                "zxgrid": False,
393            }
394        else:
395            options["axes"] = axes  # put back
396        if cmap is None:
397            cmap = "rainbow"
398    elif style in (1, 'matplotlib'):
399        bg = options.pop("bg", None)
400        if bg is None:
401            options["bg"] = "white"
402        else:
403            options["bg"] = bg
404        axes = options.pop("axes", None)
405        if axes is None:
406            options["axes"] = {
407                "xygrid": False,
408                "yzgrid": False,
409                "zxgrid": False,
410            }
411        else:
412            options["axes"] = axes  # put back
413        if cmap is None:
414            cmap = "viridis"
415    elif style in (2, 'paraview'):
416        bg = options.pop("bg", None)
417        if bg is None:
418            options["bg"] = (82, 87, 110)
419        else:
420            options["bg"] = bg
421        if cmap is None:
422            cmap = "coolwarm"
423    elif style in (3, 'meshlab'):
424        bg = options.pop("bg", None)
425        if bg is None:
426            options["bg"] = (8, 8, 16)
427            options["bg2"] = (117, 117, 234)
428        else:
429            options["bg"] = bg
430        axes = options.pop("axes", None)
431        if axes is None:
432            options["axes"] = 10
433        else:
434            options["axes"] = axes  # put back
435        if cmap is None:
436            cmap = "afmhot"
437    elif style in (4, 'bw'):
438        bg = options.pop("bg", None)
439        if bg is None:
440            options["bg"] = (217, 255, 238)
441        else:
442            options["bg"] = bg
443        axes = options.pop("axes", None)
444        if axes is None:
445            options["axes"] = {
446                "xygrid": False,
447                "yzgrid": False,
448                "zxgrid": False,
449            }
450        else:
451            options["axes"] = axes  # put back
452        if cmap is None:
453            cmap = "binary"
454
455    #################################################################
456    actors = []
457    if vedo.plotter_instance:
458        if add:
459            actors = vedo.plotter_instance.actors
460
461    if mesh and ("mesh" in mode or "color" in mode or "displace" in mode):
462
463        actor = MeshActor(u, mesh, exterior=exterior)
464
465        actor.wireframe(wire)
466        actor.scale(scaleMeshFactors)
467        if lighting:
468            actor.lighting(lighting)
469        if ttime:
470            actor.z(ttime)
471        if legend:
472            actor.legend(legend)
473        if c:
474            actor.color(c)
475        if lc:
476            actor.linecolor(lc)
477        if alpha:
478            alpha = min(alpha, 1)
479            actor.alpha(alpha * alpha)
480        if lw:
481            actor.linewidth(lw)
482            if wire and alpha:
483                lw1 = min(lw, 1)
484                actor.alpha(alpha * lw1)
485        if ps:
486            actor.pointSize(ps)
487        if shading:
488            if shading == "phong":
489                actor.phong()
490            elif shading == "flat":
491                actor.flat()
492            elif shading[0] == "g":
493                actor.gouraud()
494
495        if "displace" in mode:
496            actor.move(u)
497
498        if cmap and (actor.u_values is not None) and len(actor.u_values) and c is None:
499            if actor.u_values.ndim > 1:
500                actor.cmap(cmap, utils.mag(actor.u_values), vmin=vmin, vmax=vmax)
501            else:
502                actor.cmap(cmap, actor.u_values, vmin=vmin, vmax=vmax)
503
504        if warpYfactor:
505            scals = actor.pointdata[0]
506            if len(scals):
507                pts_act = actor.points()
508                pts_act[:, 1] = scals * warpYfactor * scaleMeshFactors[1]
509        if warpZfactor:
510            scals = actor.pointdata[0]
511            if len(scals):
512                pts_act = actor.points()
513                pts_act[:, 2] = scals * warpZfactor * scaleMeshFactors[2]
514        if warpYfactor or warpZfactor:
515            actor.points(pts_act)
516            if vmin is not None and vmax is not None:
517                actor.mapper().SetScalarRange(vmin, vmax)
518
519        if scbar and c is None:
520            if "3d" in scbar:
521                actor.add_scalarbar3d()
522            elif "h" in scbar:
523                actor.add_scalarbar(horizontal=True)
524            else:
525                actor.add_scalarbar(horizontal=False)
526
527        if len(isolns) > 0:
528            ison = isolns.pop("n", 10)
529            isocol = isolns.pop("c", "black")
530            isoalpha = isolns.pop("alpha", 1)
531            isolw = isolns.pop("lw", 1)
532
533            isos = actor.isolines(n=ison).color(isocol).lw(isolw).alpha(isoalpha)
534
535            isoz = isolns.pop("z", None)
536            if isoz is not None:  # kind of hack to make isolines visible on flat meshes
537                d = isoz
538            else:
539                d = actor.diagonal_size() / 400
540            isos.z(actor.z() + d)
541            actors.append(isos)
542
543        actors.append(actor)
544
545    #################################################################
546    if "streamline" in mode:
547        mode = mode.replace("streamline", "")
548        str_act = MeshStreamLines(u, **streamlines)
549        actors.append(str_act)
550
551    #################################################################
552    if "arrow" in mode or "line" in mode:
553        if "arrow" in mode:
554            arrs = MeshArrows(u, scale=scale)
555        else:
556            arrs = MeshLines(u, scale=scale)
557
558        if arrs:
559            if legend and "mesh" not in mode:
560                arrs.legend(legend)
561            if c:
562                arrs.color(c)
563                arrs.color(c)
564            if alpha:
565                arrs.alpha(alpha)
566            actors.append(arrs)
567
568    #################################################################
569    if "tensor" in mode:
570        pass  # todo
571
572    #################################################################
573    for ob in inputobj:
574        inputtype = str(type(ob))
575        if "vedo" in inputtype:
576            actors.append(ob)
577
578    if text:
579        # textact = Text2D(text, font=font)
580        actors.append(text)
581
582    if "at" in options and "interactive" not in options:
583        if vedo.plotter_instance:
584            N = vedo.plotter_instance.shape[0] * vedo.plotter_instance.shape[1]
585            if options["at"] == N - 1:
586                options["interactive"] = True
587
588    # if vedo.plotter_instance:
589    #     for a2 in vedo.collectable_actors:
590    #         if isinstance(a2, vtk.vtkCornerAnnotation):
591    #             if 0 in a2.rendered_at: # remove old message
592    #                 vedo.plotter_instance.remove(a2)
593    #                 break
594
595    if len(actors) == 0:
596        print('Warning: no objects to show, check mode in plot(mode="...")')
597
598    if returnActorsNoShow:
599        return actors
600
601    return show(actors, **options)
602
603
604###################################################################################
605class MeshActor(Mesh):
606    """MeshActor for dolfin support."""
607    def __init__(self, *inputobj, **options):
608        """MeshActor, a `vedo.Mesh` derived object for dolfin support."""
609
610        c = options.pop("c", None)
611        alpha = options.pop("alpha", 1)
612        exterior = options.pop("exterior", False)
613        compute_normals = options.pop("compute_normals", False)
614
615        mesh, u = _inputsort(inputobj)
616        if not mesh:
617            return
618
619        if exterior:
620            import dolfin
621
622            meshc = dolfin.BoundaryMesh(mesh, "exterior")
623        else:
624            meshc = mesh
625
626        if hasattr(mesh, "coordinates"):
627            coords = mesh.coordinates()
628        else:
629            coords = mesh.geometry.points
630
631        cells = meshc.cells()
632
633        if cells.shape[1] == 4:
634            # something wrong in this as it cannot reproduce the tet cell..
635            # from vedo.tetmesh import _buildtetugrid
636            # cells[:,[2, 0]] = cells[:,[0, 2]]
637            # cells[:,[1, 0]] = cells[:,[0, 1]]
638            # cells[:,[0, 1, 2, 3]] = cells[:,[0, 2, 1, 3]]
639            # cells[:,[0, 1, 2, 3]] = cells[:,[0, 2, 1, 3]]
640            # cells[:,[0, 1, 2, 3]] = cells[:,  [0, 1, 3, 2]]
641            # cells[:,[0, 1, 2, 3]] = cells[:,[1, 0, 2, 3]]
642            # cells[:,[0, 1, 2, 3]] = cells[:,[2, 0, 1, 3]]
643            # cells[:,[0, 1, 2, 3]] = cells[:,[2, 0, 1, 3]]
644            # print(cells[0])
645            # print(coords[cells[0]])
646            # poly = utils.geometry(_buildtetugrid(coords, cells))
647            # poly = utils.geometry(vedo.TetMesh([coords, cells]).inputdata())
648
649            poly = vtk.vtkPolyData()
650
651            source_points = vtk.vtkPoints()
652            source_points.SetData(utils.numpy2vtk(coords, dtype=np.float32))
653            poly.SetPoints(source_points)
654
655            source_polygons = vtk.vtkCellArray()
656            for f in cells:
657                # do not use vtkTetra() because it fails
658                # with dolfin faces orientation
659                ele0 = vtk.vtkTriangle()
660                ele1 = vtk.vtkTriangle()
661                ele2 = vtk.vtkTriangle()
662                ele3 = vtk.vtkTriangle()
663
664                f0, f1, f2, f3 = f
665                pid0 = ele0.GetPointIds()
666                pid1 = ele1.GetPointIds()
667                pid2 = ele2.GetPointIds()
668                pid3 = ele3.GetPointIds()
669
670                pid0.SetId(0, f0)
671                pid0.SetId(1, f1)
672                pid0.SetId(2, f2)
673
674                pid1.SetId(0, f0)
675                pid1.SetId(1, f1)
676                pid1.SetId(2, f3)
677
678                pid2.SetId(0, f1)
679                pid2.SetId(1, f2)
680                pid2.SetId(2, f3)
681
682                pid3.SetId(0, f2)
683                pid3.SetId(1, f3)
684                pid3.SetId(2, f0)
685
686                source_polygons.InsertNextCell(ele0)
687                source_polygons.InsertNextCell(ele1)
688                source_polygons.InsertNextCell(ele2)
689                source_polygons.InsertNextCell(ele3)
690
691            poly.SetPolys(source_polygons)
692
693        else:
694            poly = utils.buildPolyData(coords, cells)
695
696        Mesh.__init__(
697            self,
698            poly,
699            c=c,
700            alpha=alpha,
701        )
702        if compute_normals:
703            self.compute_normals()
704
705        self.mesh = mesh  # holds a dolfin Mesh obj
706        self.u = u  # holds a dolfin function_data
707        # holds the actual values of u on the mesh
708        self.u_values = _compute_uvalues(u, mesh)
709
710    def move(self, u=None, deltas=None):
711        """Move mesh according to solution `u` or from calculated vertex displacements `deltas`."""
712        if u is None:
713            u = self.u
714        if deltas is None:
715            if self.u_values is not None:
716                deltas = self.u_values
717            else:
718                deltas = _compute_uvalues(u, self.mesh)
719                self.u_values = deltas
720
721        if hasattr(self.mesh, "coordinates"):
722            coords = self.mesh.coordinates()
723        else:
724            coords = self.mesh.geometry.points
725
726        if coords.shape != deltas.shape:
727            vedo.logger.error(
728                f"Try to move mesh with wrong solution type shape {coords.shape} vs {deltas.shape}"
729            )
730            vedo.logger.error("Mesh is not moved. Try mode='color' in plot().")
731            return
732
733        movedpts = coords + deltas
734        if movedpts.shape[1] == 2:  # 2d
735            movedpts = np.c_[movedpts, np.zeros(movedpts.shape[0])]
736        self.polydata(False).GetPoints().SetData(utils.numpy2vtk(movedpts, dtype=np.float32))
737        self.polydata(False).GetPoints().Modified()
738
739
740def MeshPoints(*inputobj, **options):
741    """Build a point object of type `Mesh` for a list of points."""
742    r = options.pop("r", 5)
743    c = options.pop("c", "gray")
744    alpha = options.pop("alpha", 1)
745
746    mesh, u = _inputsort(inputobj)
747    if not mesh:
748        return None
749
750    if hasattr(mesh, "coordinates"):
751        plist = mesh.coordinates()
752    else:
753        plist = mesh.geometry.points
754
755    u_values = _compute_uvalues(u, mesh)
756
757    if len(plist[0]) == 2:  # coords are 2d.. not good..
758        plist = np.insert(plist, 2, 0, axis=1)  # make it 3d
759    if len(plist[0]) == 1:  # coords are 1d.. not good..
760        plist = np.insert(plist, 1, 0, axis=1)  # make it 3d
761        plist = np.insert(plist, 2, 0, axis=1)
762
763    actor = shapes.Points(plist, r=r, c=c, alpha=alpha)
764
765    actor.mesh = mesh
766    if u:
767        actor.u = u
768        if len(u_values.shape) == 2:
769            if u_values.shape[1] in [2, 3]:  # u_values is 2D or 3D
770                actor.u_values = u_values
771                dispsizes = utils.mag(u_values)
772        else:  # u_values is 1D
773            dispsizes = u_values
774        actor.pointdata["u_values"] = dispsizes
775        actor.pointdata.select("u_values")
776    return actor
777
778
779def MeshLines(*inputobj, **options):
780    """
781    Build the line segments between two lists of points `start_points` and `end_points`.
782    `start_points` can be also passed in the form `[[point1, point2], ...]`.
783
784    A dolfin `Mesh` that was deformed/modified by a function can be
785    passed together as inputs.
786
787    Use `scale` to apply a rescaling factor to the length
788    """
789    scale = options.pop("scale", 1)
790    lw = options.pop("lw", 1)
791    c = options.pop("c", "grey")
792    alpha = options.pop("alpha", 1)
793
794    mesh, u = _inputsort(inputobj)
795    if not mesh:
796        return None
797
798    if hasattr(mesh, "coordinates"):
799        start_points = mesh.coordinates()
800    else:
801        start_points = mesh.geometry.points
802
803    u_values = _compute_uvalues(u, mesh)
804    if not utils.is_sequence(u_values[0]):
805        vedo.logger.error("cannot show Lines for 1D scalar values")
806        raise RuntimeError()
807
808    end_points = start_points + u_values
809    if u_values.shape[1] == 2:  # u_values is 2D
810        u_values = np.insert(u_values, 2, 0, axis=1)  # make it 3d
811        start_points = np.insert(start_points, 2, 0, axis=1)  # make it 3d
812        end_points = np.insert(end_points, 2, 0, axis=1)  # make it 3d
813
814    actor = shapes.Lines(start_points, end_points, scale=scale, lw=lw, c=c, alpha=alpha)
815
816    actor.mesh = mesh
817    actor.u = u
818    actor.u_values = u_values
819    return actor
820
821
822def MeshArrows(*inputobj, **options):
823    """Build arrows representing displacements."""
824    s = options.pop("s", None)
825    c = options.pop("c", "gray")
826    scale = options.pop("scale", 1)
827    alpha = options.pop("alpha", 1)
828    res = options.pop("res", 12)
829
830    mesh, u = _inputsort(inputobj)
831    if not mesh:
832        return None
833
834    if hasattr(mesh, "coordinates"):
835        start_points = mesh.coordinates()
836    else:
837        start_points = mesh.geometry.points
838
839    u_values = _compute_uvalues(u, mesh)
840    if not utils.is_sequence(u_values[0]):
841        vedo.logger.error("cannot show Arrows for 1D scalar values")
842        raise RuntimeError()
843
844    end_points = start_points + u_values * scale
845    if u_values.shape[1] == 2:  # u_values is 2D
846        u_values = np.insert(u_values, 2, 0, axis=1)  # make it 3d
847        start_points = np.insert(start_points, 2, 0, axis=1)  # make it 3d
848        end_points = np.insert(end_points, 2, 0, axis=1)  # make it 3d
849
850    actor = shapes.Arrows(start_points, end_points, s=s, alpha=alpha, res=res)
851    actor.color(c)
852    actor.mesh = mesh
853    actor.u = u
854    actor.u_values = u_values
855    return actor
856
857
858def MeshStreamLines(*inputobj, **options):
859    """Build a streamplot."""
860    from vedo.shapes import StreamLines
861
862    print("Building streamlines...")
863
864    tol = options.pop("tol", 0.02)
865    lw = options.pop("lw", 2)
866    direction = options.pop("direction", "forward")
867    max_propagation = options.pop("max_propagation", None)
868    scalar_range = options.pop("scalar_range", None)
869    probes = options.pop("probes", None)
870
871    tubes = options.pop("tubes", {})  # todo
872    maxRadiusFactor = options.pop("maxRadiusFactor", 1)
873    varyRadius = options.pop("varyRadius", 1)
874
875    mesh, u = _inputsort(inputobj)
876    if not mesh:
877        return None
878
879    u_values = _compute_uvalues(u, mesh)
880    if not utils.is_sequence(u_values[0]):
881        vedo.logger.error("cannot show Arrows for 1D scalar values")
882        raise RuntimeError()
883    if u_values.shape[1] == 2:  # u_values is 2D
884        u_values = np.insert(u_values, 2, 0, axis=1)  # make it 3d
885
886    meshact = MeshActor(u)
887    meshact.pointdata["u_values"] = u_values
888    meshact.pointdata.select("u_values")
889
890    if utils.is_sequence(probes):
891        pass  # it's already it
892    elif tol:
893        print("decimating mesh points to use them as seeds...")
894        probes = meshact.clone().subsample(tol).points()
895    else:
896        probes = meshact.points()
897    if len(probes) > 500:
898        printc("Probing domain with n =", len(probes), "points")
899        printc(" ..this may take time (or choose a larger tol value)")
900
901    if lw:
902        tubes = {}
903    else:
904        tubes["varyRadius"] = varyRadius
905        tubes["maxRadiusFactor"] = maxRadiusFactor
906
907    str_lns = StreamLines(
908        meshact,
909        probes,
910        direction=direction,
911        max_propagation=max_propagation,
912        tubes=tubes,
913        scalar_range=scalar_range,
914        active_vectors="u_values",
915    )
916
917    if lw:
918        str_lns.lw(lw)
919
920    return str_lns