# Visualize non-Lagrange vector fields

**URL:** <https://fenicsproject.discourse.group/t/visualize-non-lagrange-vector-fields/15663>\
**Category:** dolfinx\
**Created:** [August 28, 2024, 12:20pm UTC](https://fenicsproject.discourse.group/t/visualize-non-lagrange-vector-fields/15663 "2024-08-28T12:20:14Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![andrmykl](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/andrmykl/32/7481_2.png) [@andrmykl](https://fenicsproject.discourse.group/u/andrmykl)\
**Post date:** [August 28, 2024, 12:20pm UTC](https://fenicsproject.discourse.group/t/visualize-non-lagrange-vector-fields/15663/1 "2024-08-28T12:20:14Z")

</div>

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?

 ![Screenshot from 2024-08-28 14-17-37](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/8/820011c5effc9105aa32bd290465f8db67c7d267.jpeg)

```auto
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()

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 28, 2024, 12:43pm UTC](https://fenicsproject.discourse.group/t/visualize-non-lagrange-vector-fields/15663/2 "2024-08-28T12:43:30Z")

</div>

> [@andrmykl](#):
>
> `P = basix.ufl.element(basix.ElementFamily.P,basix.CellType.triangle, shape = (3,), degree = 2, discontinuous=False)`

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

---

<div class="post-metadata">

**Author:** ![andrmykl](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/andrmykl/32/7481_2.png) [@andrmykl](https://fenicsproject.discourse.group/u/andrmykl)\
**Post date:** [August 28, 2024, 1:06pm UTC](https://fenicsproject.discourse.group/t/visualize-non-lagrange-vector-fields/15663/3 "2024-08-28T13:06:08Z")

</div>

These are the plots when I use the discontinuous version,

 ![Screenshot from 2024-08-28 15-04-48](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/9/9628ff072c36dd321ee4adfb7d898b1bb01cea25.jpeg)

which does not seen to help much.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 28, 2024, 1:23pm UTC](https://fenicsproject.discourse.group/t/visualize-non-lagrange-vector-fields/15663/4 "2024-08-28T13:23:37Z")

</div>

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:

```python
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

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

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 28, 2024, 1:38pm UTC](https://fenicsproject.discourse.group/t/visualize-non-lagrange-vector-fields/15663/5 "2024-08-28T13:38:05Z")

</div>

Issue submitted at: [[BUG]: Interpolation on manifold broken · Issue #3372 · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/issues/3372)

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 28, 2024, 3:11pm UTC](https://fenicsproject.discourse.group/t/visualize-non-lagrange-vector-fields/15663/6 "2024-08-28T15:11:08Z")

</div>

Issue partially resolved with:

> <https://github.com/FEniCS/dolfinx/pull/3373>
>
> Use element value size to create reference buffer.
> At least partially resolves …#3372

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:

> **[DefElement](https://defelement.com/elements/examples/triangle-nedelec1-lagrange-1.html)**
>
> An encyclopedia of finite element definitions

Visualizing the interpolation into P of the pure function yields

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/41facc978c83a189400288c31d40dc59fe163f57.png)  
and in the x-y plane looks like:  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/d97f56e9a654adea991ce712777238ea48d929c9.png)  
However, when you interpolate `f` into N1E, you get  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/b/bc465f5c7c2d504320e0b50f48c6e9ccf63b748e.png)  
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.

```python
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))

```