vedo.dolfin

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