Visualize non-Lagrange vector fields

andrmykl · August 28, 2024, 12:20pm

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

dokken · August 28, 2024, 12:43pm

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

andrmykl · August 28, 2024, 1:06pm

These are the plots when I use the discontinuous version,

which does not seen to help much.

dokken · August 28, 2024, 1:23pm

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

dokken · August 28, 2024, 1:38pm

Issue submitted at: [BUG]: Interpolation on manifold broken · Issue #3372 · FEniCS/dolfinx · GitHub

dokken · August 28, 2024, 3:11pm

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

Topic		Replies	Views
3D vector field: plot 2D slice of 3D mesh dolfinx mesh	5	228	April 12, 2024
Dolfinx and pygmsh : error when interpolating dolfinx	1	771	August 21, 2020
Identifying cells and their vertices mesh	7	1605	November 4, 2022
Plotting reaction forces on a mesh General	1	54	June 25, 2024
Submesh from parallel mesh dolfinx	9	40	September 5, 2024

Visualize non-Lagrange vector fields

Related Topics