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