# Visualize non-Lagrange vector fields

I want to plot vector fields belonging to spaces other than Lagrange spaces. As far as I know, the way to do this is to interpolate to a Lagrange vector field, and plot the Lagrange vector field. I have tried to do this, but it results in some weird plots. What is going on? And is this the right way to plot vector fields that are not Lagrange?

``````import numpy as np
import gmsh
from dolfinx.io import gmshio
from dolfinx import plot
import basix
from mpi4py import MPI
import pyvista as pv
import dolfinx.mesh as mesh
import dolfinx.fem as fem

approximate_simplex_diameter = 0.25
deviation = 0.001

gmsh.initialize()
gmsh.option.setNumber("General.Terminal",0) # Turn off terminal output for gmsh

model = gmsh.model()
model.add("Sphere") # Adds a new model with the name "name"
model.setCurrent("Sphere") # Sets the model we are currently working on to the model called "name"
sphere = model.occ.addSphere(0,0,0,1) # Adds a sphere with radius 1 and center at (0,0,0) to the model
centre_of_earth = model.occ.addPoint(0,0,0) # Adds a point at the center of the sphere
point_close_to_point2 = model.occ.addPoint(np.cos(-np.pi/2 -approximate_simplex_diameter),0,np.sin(-np.pi/2 -approximate_simplex_diameter)) # Adds a point close to the south pole
equator = model.occ.addCircleArc(1,centre_of_earth,point_close_to_point2)

gmsh.model.occ.synchronize()
model.add_physical_group(dim=2, tags=[sphere]) # Adds a physical group to the model. The group is a 2D group with the tag "sphere"
gmsh.model.mesh.embed(1, [equator],2,1) # Embeds the 1D entities in the 2D entities

gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.option.setNumber("Mesh.ElementOrder", 2)
#gmsh.option.setNumber("Mesh.HighOrderOptimize", 2)
gmsh.option.setNumber("Mesh.Smoothing", 100)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", approximate_simplex_diameter - deviation)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", approximate_simplex_diameter + deviation)
gmsh.model.mesh.generate(2)

msh, cell_tags,facet_tags= gmshio.model_to_mesh(model,MPI.COMM_WORLD, gdim = 3, rank = 0)

N1E = basix.ufl.element(basix.ElementFamily.N1E,basix.CellType.triangle,degree = 1)
BDM = basix.ufl.element(basix.ElementFamily.BDM,basix.CellType.triangle, degree = 1)
P = basix.ufl.element(basix.ElementFamily.P,basix.CellType.triangle, shape = (3,), degree = 2, discontinuous=False)

N1E = fem.functionspace(msh, N1E)
BDM = fem.functionspace(msh, BDM)
P = fem.functionspace(msh, P)

n1e = fem.Function(N1E)
bdm = fem.Function(BDM)
p = fem.Function(P)

n1e.interpolate(lambda x : np.array([-x[2], np.zeros_like(x[1]), x[0]]))
bdm.interpolate(lambda x : np.array([ -x[2], np.zeros_like(x[1]), x[0]]))
p.interpolate(lambda x : np.array([ -x[2], np.zeros_like(x[1]), x[0]]))

cells, cell_types, x = plot.vtk_mesh(P)
grid = pv.UnstructuredGrid(cells, cell_types, x)

plotter = pv.Plotter(shape=(1, 3))
plotter.subplot(0, 0)
plotter.add_mesh(grid,style = "surface", color = "beige")
grid.point_data["v"] = p.x.array.reshape(-1, msh.geometry.dim)
plotter.add_text("Interpolate directly", font_size=12, position="upper_edge")
glyphs = grid.glyph(orient="v", factor=0.08)
plotter.add_mesh(glyphs, show_scalar_bar=False)

p.interpolate(n1e)
plotter.subplot(0, 1)
plotter.add_mesh(grid,style = "surface" ,color = "beige")
grid.point_data["v"] = p.x.array.reshape(-1, msh.geometry.dim)
plotter.add_text("Interpolate via N1E", font_size=12, position="upper_edge")
glyphs = grid.glyph(orient="v", factor=0.08)
plotter.add_mesh(glyphs, show_scalar_bar=False)

p.interpolate(bdm)
plotter.subplot(0, 2)
plotter.add_mesh(grid,style = "surface" ,color = "beige")
grid.point_data["v"] = p.x.array.reshape(-1, msh.geometry.dim)
plotter.add_text("Interpolate via BDM", font_size=12, position="upper_edge")
glyphs = grid.glyph(orient="v", factor=0.08)
plotter.add_mesh(glyphs, show_scalar_bar=False)

plotter.show()
``````

You need to use a discontinuous Lagrange element if want to represent functions from N1E or BDM, as they only have partial inter-element continuity.
Change the `discontinuous`-flag to `True`

Edit
I am also not certain about what polynomial spaces that are embedded in BDM or N1E or a manifold

These are the plots when I use the discontinuous version,

which does not seen to help much.

One thing that seems a bit odd in your example:
You are using a second order mesh for the manifold, but a first order function space for N1E and BDM.

However, using a second order space doesn’t reduce the error:

``````import numpy as np
import gmsh
from dolfinx.io import gmshio
import basix
from mpi4py import MPI
import dolfinx.fem as fem

approximate_simplex_diameter = 0.25
deviation = 0.001

gmsh.initialize()
gmsh.option.setNumber("General.Terminal",0) # Turn off terminal output for gmsh

model = gmsh.model()
model.add("Sphere") # Adds a new model with the name "name"
model.setCurrent("Sphere") # Sets the model we are currently working on to the model called "name"
sphere = model.occ.addSphere(0,0,0,1) # Adds a sphere with radius 1 and center at (0,0,0) to the model
centre_of_earth = model.occ.addPoint(0,0,0) # Adds a point at the center of the sphere
point_close_to_point2 = model.occ.addPoint(np.cos(-np.pi/2 -approximate_simplex_diameter),0,np.sin(-np.pi/2 -approximate_simplex_diameter)) # Adds a point close to the south pole
equator = model.occ.addCircleArc(1,centre_of_earth,point_close_to_point2)

gmsh.model.occ.synchronize()
model.add_physical_group(dim=2, tags=[sphere]) # Adds a physical group to the model. The group is a 2D group with the tag "sphere"
gmsh.model.mesh.embed(1, [equator],2,1) # Embeds the 1D entities in the 2D entities

gmsh.option.setNumber("Mesh.Algorithm", 6)
degree = 2
gmsh.option.setNumber("Mesh.ElementOrder", degree)
#gmsh.option.setNumber("Mesh.HighOrderOptimize", 2)
gmsh.option.setNumber("Mesh.Smoothing", 100)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", approximate_simplex_diameter - deviation)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", approximate_simplex_diameter + deviation)
gmsh.model.mesh.generate(2)

msh, cell_tags,facet_tags= gmshio.model_to_mesh(model,MPI.COMM_WORLD, gdim = 3, rank = 0)

N1E = basix.ufl.element(basix.ElementFamily.N1E,basix.CellType.triangle,degree = degree)
BDM = basix.ufl.element(basix.ElementFamily.BDM,basix.CellType.triangle, degree = degree)
P = basix.ufl.element(basix.ElementFamily.P,basix.CellType.triangle, shape = (3,), degree = degree, discontinuous=True)

N1E = fem.functionspace(msh, N1E)
BDM = fem.functionspace(msh, BDM)
P = fem.functionspace(msh, P)

n1e = fem.Function(N1E)
bdm = fem.Function(BDM)
p = fem.Function(P)

def f(x):
return 0*x[0], 0*x[1], 1*x[2]

n1e.interpolate(f)
bdm.interpolate(f)
p.interpolate(f)

import ufl
def error(ph):
x_ = ufl.SpatialCoordinate(msh)
f_ex = ufl.as_vector(f(x_))
diff = ufl.inner(ph - f_ex, ph-f_ex) * ufl.dx
compiled = fem.form(diff)
local_err = fem.assemble_scalar(compiled)
return np.sqrt(MPI.COMM_WORLD.allreduce(local_err, op=MPI.SUM))

print(f"{error(p)=}, {error(bdm)=}, {error(n1e)=}")
``````

yielding

``````rror(p)=2.654838473666969e-16, error(bdm)=2.043703290434796, error(n1e)=36.528579395455964

``````

Issue partially resolved with:

However, consider the following MWE.
I have a single skewed triangle, which I interpolate (-y, x) into, which is one of the basis functions for a flat N1E element:

Visualizing the interpolation into P of the pure function yields

and in the x-y plane looks like:

However, when you interpolate `f` into N1E, you get

which aligns perfectly with the surface tangents of the cell, which gets you to the point:

• For manifold geometries, standard polynomials might not be representable in the function space. Maybe @jpdean or @mscroggs could sanity-check this logic.
``````from mpi4py import MPI

import numpy as np
import basix.ufl
import ufl
import dolfinx as dfx
def error(ph):
x_ = ufl.SpatialCoordinate(msh)
f_ex = ufl.as_vector(f(x_))
diff = ufl.inner(ph - f_ex, ph-f_ex) * ufl.dx
compiled = dfx.fem.form(diff)
local_err = dfx.fem.assemble_scalar(compiled)
return np.sqrt(MPI.COMM_WORLD.allreduce(local_err, op=MPI.SUM))

gdim = 3

vertices = np.array(
[(0.1, 0.2, 1), (0.2, 1.1, 3), (1.3, 1.4, 3)],
dtype=dfx.default_real_type,
)
vertices = vertices[:, :gdim]
cells = [(0, 1, 2)]
domain = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(gdim,), dtype=dfx.default_real_type))
msh = dfx.mesh.create_mesh(MPI.COMM_WORLD, cells, vertices, domain)

N1E = basix.ufl.element(basix.ElementFamily.N1E,basix.CellType.triangle,degree=1)

N1E = dfx.fem.functionspace(msh, N1E)

n1e = dfx.fem.Function(N1E)

def f(x):
if gdim == 2:
return -x[1], x[0]
elif gdim == 3:
return -x[1], x[0], 0*x[2]

n1e.interpolate(f)
print(f"N1e {gdim=} {error(n1e)=}")
degree=2
P = basix.ufl.element(basix.ElementFamily.P,basix.CellType.triangle, shape = (gdim,), degree = degree, discontinuous=True)

P = dfx.fem.functionspace(msh, P)
p = dfx.fem.Function(P)

with dfx.io.VTXWriter(msh.comm, "uh.bp", [p]) as bp:
p.interpolate(f)
bp.write(0.0)

#print(p.x.array)
p.interpolate(n1e)
bp.write(1.0)
print(error(p))
``````