3D vector field: plot 2D slice of 3D mesh

dolfinx
1

Hi all,

I’m solving a nonlinear problem in a 3D mesh, the result of which is a 3D vector field:

import numpy as np
import ufl
import gmsh
from dolfinx import fem, nls
from dolfinx.io import gmshio
from petsc4py.PETSc import ScalarType
from mpi4py import MPI

model = gmsh.model; geom = gmsh.model.geo
mesh_comm = MPI.COMM_WORLD

def generateSphere3D(RR):
    """Create a Gmsh model of a sphere."""
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    model.add("model_acorn")
    sphere = model.occ.addSphere(0, 0, 0, RR, tag=1)
    model.occ.synchronize()
    model.add_physical_group(dim=3, tags=[sphere])    
    geom.synchronize()
    model.mesh.generate(3)
    return gmshio.model_to_mesh(model, mesh_comm, 0, gdim=3)

domain, _, _ = generateSphere3D(1)

#### define the space of funtions and the functions that we will use
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, 3)
V = fem.FunctionSpace(domain, element) # space of vector functions with 3 components
m = fem.Function(V, name='m') # 3-dimensional vector OP
dm = ufl.TrialFunction(V) # variation of the OP
phi = ufl.TestFunction(V) # test function

#### derive the weak formulation directly from the energy functional
def ET(u):
    # just a simple example
    a = -1; c = 1
    return (a*ufl.dot(u,u)/2. + c*ufl.dot(u,u)**2/4.)*ufl.dx    
F = ET(m)
dF = ufl.derivative(F, m, phi) # variation of the energy
JF = ufl.derivative(dF, m, dm) # Jacobian of the energy

#### define and solve the problem
problem = fem.petsc.NonlinearProblem(dF, m, J=JF)

def solveProblem(prob):
    solver = nls.petsc.NewtonSolver(mesh_comm, prob)
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-6 
    solver.max_it = 100
    solver.report = True
    minit = [1.,0,0] # initial guess for the unknown function
    minit = ufl.as_vector([fem.Constant(domain, ScalarType(minit[0])),
                           fem.Constant(domain, ScalarType(minit[1])),
                           fem.Constant(domain, ScalarType(minit[2]))])
    expr = fem.Expression(minit, V.element.interpolation_points())
    m.interpolate(expr)
    solver.nonzero_initial_guess = True
    n, converged = solver.solve(m)
    assert(converged)

solveProblem(problem)

Now, for visualisation purposes, I want to plot the 3D-vector solution for the y=0 slice (2D) of the original system. Building on this link, I tried this:

def generateCircle2D(RR, ms=0.1):
    """Create a Gmsh model of a circle."""
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    model.add("y=0 projection")
    p = [
        geom.addPoint(0, 0, RR, ms),
        geom.addPoint(0, 0, 0, ms),
        geom.addPoint(RR, 0, 0, ms),
        geom.addPoint(0, 0, -RR, ms),
        geom.addPoint(-RR, 0, 0, ms)
    ]
    l = [
        geom.addCircleArc(p[0], p[1], p[2]),
        geom.addCircleArc(p[2], p[1], p[3]),
        geom.addCircleArc(p[3], p[1], p[4]),
        geom.addCircleArc(p[4], p[1], p[0])
    ]
    cl = geom.addCurveLoop([l[j] for j in range(len(l))])
    s = geom.addPlaneSurface([cl]) 
    model.occ.synchronize()
    model.add_physical_group(2, [s])
    geom.synchronize()
    model.mesh.generate(2)
    return gmshio.model_to_mesh(model, mesh_comm, 0, gdim=3)

domain2D, _, _ = generateCircle2D(1)
element2D = ufl.VectorElement("CG", domain2D.ufl_cell(), 1, 3)
V2D = fem.FunctionSpace(domain2D, element2D)

class m_slice(fem.Expression):
    def eval(self, value, x):
        value[0] = m(x[0], 0, x[2]) # slice of 3D system: y = 0 plane
    def value_shape(self):
        return ()

m_sl = fem.Function.interpolate(m_slice, V2D)

However, I get the error message:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [5], in <cell line: 35>()
     32     def value_shape(self):
     33         return ()
---> 35 m_sl = fem.Function.interpolate(m_slice, V2D)

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py:333, in Function.interpolate(self, u, cells)
    330     self._cpp_object.interpolate(expr._cpp_object, cells)
    332 if cells is None:
--> 333     mesh = self.function_space.mesh
    334     map = mesh.topology.index_map(mesh.topology.dim)
    335     cells = np.arange(map.size_local + map.num_ghosts, dtype=np.int32)

AttributeError: type object 'm_slice' has no attribute 'function_space'

The idea is then to generate a stream plot of the projected solution:

from scipy.interpolate import griddata
def plotStreamlines(mm, res=100):
    """Plots the OP streamlines and magnitude in Cartesian coordinates."""
    
    m_array = mm.x.array 
    m_array = [ [ m_array[i], m_array[i+1], m_array[i+2] ] for i in range(0, len(m_array), 3) ]
    m_array = np.array(m_array).T
    x_array = V2D.tabulate_dof_coordinates() 
    x_array = np.array(x_array).T
    x_grid, z_grid = np.meshgrid(np.linspace(x_array[0].min(),x_array[0].max(),res),
                                 np.linspace(x_array[2].min(),x_array[2].max(),res))   
    x_vector_interp = griddata((x_array[0], x_array[2]), m_array[0], (x_grid, z_grid), method='cubic')
    z_vector_interp = griddata((x_array[0], x_array[2]), m_array[2], (x_grid, z_grid), method='cubic')
    
    plt.streamplot(x_grid, z_grid, x_vector_interp, z_vector_interp)
    plt.show()

plotStreamlines(m_sl)

Could anyone please help me? Thanks in advance!

2

I think it is best to use dedicated visualization tools for that, such as Paraview for instance

2 Likes
3

Agreed, it will be quite a task to cut a slice of a mesh from a finite element perspective (with arbitrary elements), with very little gain. I would use Paraview or Pyvista for this.

4

Dear @bleyerj and @dokken,
thanks for your replies. I could make a 3D plot of the solution with this function:

def plotPyvista3D():
    pyvista.set_jupyter_backend("pythreejs")
    topology, cell_types, geometry = create_vtk_mesh(V)
    values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
    values[:, :len(m)] = m.x.array.real.reshape((geometry.shape[0], len(m)))

    # Create a point cloud of glyphs
    function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
    function_grid["u"] = values
    glyphs = function_grid.glyph(orient="u", factor=0.2)

    # Create a pyvista-grid for the mesh
    grid = pyvista.UnstructuredGrid(*create_vtk_mesh(domain, domain.topology.dim))

    # Create plotter
    plotter = pyvista.Plotter()
    plotter.add_mesh(grid, style="wireframe", color="k")
    plotter.add_mesh(glyphs)
    plotter.view_xy()
    if not pyvista.OFF_SCREEN:
        plotter.show()
    else:
        pyvista.start_xvfb()
        fig_as_array = plotter.screenshot("glyphs.png")

Could you please advise on how to plot the y=0 slice? Doing domain.slice(normal=[0,1,0]) returns AttributeError: 'Mesh' object has no attribute 'slice'.
Also, is there any way to force all the vectors to have the same length, with their norm encoded only in the color?
Thanks!

5

You are now mixing the dolfinx mesh (domain) with the pyvista.UnstructuredGrid (grid).

7

Thanks. Now I have this function, but I don’t know how to plot only the vectors that correspond to the slice:

def plotPyvistaSlice():
    pyvista.set_jupyter_backend("pythreejs")
    topology, cell_types, geometry = create_vtk_mesh(V)
    
    values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
    values[:, :len(m)] = m.x.array.real.reshape((geometry.shape[0], len(m)))
    values = [vec/np.linalg.norm(vec) for vec in values]
    
    # Create a point cloud of glyphs
    function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
    function_grid["u"] = values
    glyphs = function_grid.glyph(orient="u", factor=0.2)

    # Create a pyvista-grid for the mesh
    grid = pyvista.UnstructuredGrid(*create_vtk_mesh(domain, domain.topology.dim))
    single_slice = grid.slice(normal=[0, 1, 0])

    # Create plotter
    plotter = pyvista.Plotter()
    plotter.add_mesh(grid.outline(), color="k")
    plotter.add_mesh(single_slice, style="wireframe", color="k")
    plotter.add_mesh(glyphs, color="k")
    plotter.view_xz()
    if not pyvista.OFF_SCREEN:
        plotter.show()
    else:
        pyvista.start_xvfb()
        fig_as_array = plotter.screenshot("glyphs.png")

Anyway, my goal is actually more complex: I want to plot the normalised vectors on the slice and at the same time plot the norm with a colormap. Something like this, which I could manage using matplotlib (hence my first attempt) from a 2D solution:
image

Could you maybe help me with the code to make something similar with pyvista?

UPDATE:
I have been able to plot the vectors in the slice by doing:

def plotPyvistaSlice(norm=False):
    pyvista.set_jupyter_backend("pythreejs")
    topology, cell_types, geometry = create_vtk_mesh(V)

    # Extract the solution values
    values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
    values[:, :len(m)] = m.x.array.real.reshape((geometry.shape[0], len(m)))
    if norm: values = [vec/np.linalg.norm(vec) for vec in values]

    # Extract mesh coordinates
    coordinates = V.tabulate_dof_coordinates()

    # Find indices of points lying in the y=0 slice
    y_coords = coordinates[:, 1]
    y0_indices = np.where(np.less_equal(np.abs(y_coords), 0.02))[0]
    # print(y0_indices)

    # Extract points and corresponding values lying in the y=0 slice
    slice_points = [coordinates[ind] for ind in y0_indices]
    slice_values = [values[ind] for ind in y0_indices]

    # Create a pyvista point cloud of glyphs for the slice
    points = pyvista.PolyData(slice_points)
    points["u"] = slice_values
    glyphs = points.glyph(orient="u", factor=0.2)
    
    # grid = pyvista.UnstructuredGrid(*create_vtk_mesh(domain, domain.topology.dim))
    # single_slice = grid.slice(normal=[0, 1, 0])

    # Create plotter
    plotter = pyvista.Plotter()
    plotter.add_mesh(glyphs)
    # plotter.add_mesh(single_slice, color="grey")
    plotter.view_xz()
    plotter.set_background('white')

    plotter.show()