Split gradient of solution into vector components

General
,
1

According to this tutorial, the general idea of how to post-process a solution is to leverage dolfinx.fem.Expression(). What I want to do with my solution is find the gradient of my solution, then split up that gradient into its orthogonal vector components pointed in the \hat{\textbf{i}} and \hat{\textbf{j}} directions, and lastly use them to make a vector plot (I’ll know if what I have is correct if the vectors point normal to the gradient).

So, I need to know

  1. What expression will represent the gradient?
  2. How to split up the gradient into its vectoral components?

For 1, I tried using ufl.grad()

W = df.fem.functionspace(domain, ("DQ", 1))
expr = df.fem.Expression(ufl.grad(uh),W.element.interpolation_points())
u_grad = df.fem.Function(W)
u_grad.interpolate(expr)

And this did not work.

Here’s the setup and solution I’m working with:

# Create domain
domain = df.mesh.create_unit_square(MPI.COMM_WORLD, 61, 61, df.mesh.CellType.quadrilateral)

# Define (continuous) functionspace
V = df.fem.functionspace(domain, ("Lagrange", 1))

# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = df.mesh.exterior_facet_indices(domain.topology)

# Locate boundaries
boundary_dofs = df.fem.locate_dofs_topological(V, fdim, boundary_facets)
# Apply boundary conditions 
bc = df.fem.dirichletbc(df.default_scalar_type(0), boundary_dofs, V)

# Define Trial and Test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Define source
f = df.fem.Constant(domain, df.default_scalar_type(4))

# Variational form
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
2

If your solution uh is a scalar function living in CG1, its gradient is a vectorial function living in DG0. You should define such a space and then you can access its components through function.sub(0 or 1).collapse() like the following MWE:

import dolfinx as df
import ufl
from mpi4py import MPI
from dolfinx.fem.petsc import LinearProblem

# Create domain
domain = df.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, df.mesh.CellType.quadrilateral)

# Define (continuous) functionspace
V = df.fem.functionspace(domain, ("Lagrange", 1))

# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = df.mesh.exterior_facet_indices(domain.topology)

# Locate boundaries
boundary_dofs = df.fem.locate_dofs_topological(V, fdim, boundary_facets)
# Apply boundary conditions 
bc = df.fem.dirichletbc(df.default_scalar_type(0), boundary_dofs, V)

# Define Trial and Test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Define source
f = df.fem.Constant(domain, df.default_scalar_type(4))

# Variational form
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

W = df.fem.functionspace(domain, ("DG", 0, (domain.topology.dim, )))
expr = df.fem.Expression(ufl.grad(uh),W.element.interpolation_points())
u_grad = df.fem.Function(W)
u_grad.interpolate(expr)

u_grad_x = u_grad.sub(0).collapse()  # x component
print(u_grad_x.x.array)
u_grad_y = u_grad.sub(1).collapse()  # y component
print(u_grad_y.x.array)
2 Likes
3

Looks like it worked. Now I just need to verify the results. Would you by chance know how to plot the vectors?

4

You can export the data to visualize it in ParaView (there are many examples available on the forum), or plot it manually using the following MWE

import dolfinx as df
import ufl
from mpi4py import MPI
from dolfinx.fem.petsc import LinearProblem
import numpy as np
import matplotlib.pyplot as plt
# Create domain
domain = df.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, df.mesh.CellType.quadrilateral)

# Define (continuous) functionspace
V = df.fem.functionspace(domain, ("Lagrange", 1))

# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = df.mesh.exterior_facet_indices(domain.topology)

# Locate boundaries
boundary_dofs = df.fem.locate_dofs_topological(V, fdim, boundary_facets)
# Apply boundary conditions 
bc = df.fem.dirichletbc(df.default_scalar_type(0), boundary_dofs, V)

# Define Trial and Test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Define source
f = df.fem.Constant(domain, df.default_scalar_type(4))

# Variational form
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

W = df.fem.functionspace(domain, ("DG", 0, (domain.topology.dim, )))
expr = df.fem.Expression(ufl.grad(uh),W.element.interpolation_points())
u_grad = df.fem.Function(W)
u_grad.interpolate(expr)

u_grad_x = u_grad.sub(0).collapse().x.array  # x component
u_grad_y = u_grad.sub(1).collapse().x.array  # y component

position = W.tabulate_dof_coordinates()
x_coords = position[:, 0]
y_coords = position[:, 1]


magnitude = np.sqrt(u_grad_x**2 + u_grad_y**2)
plt.figure(figsize=(10, 8))
plt.quiver(x_coords, y_coords, u_grad_x, u_grad_y, magnitude, 
           alpha=0.8, scale=10, cmap='plasma')
plt.colorbar(label='Gradient magnitude')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.axis('equal')
plt.show()

If you want to visualize any field, you could also consider using the user-friendly pyvista-for-dolfinx package
https://gitlab.com/fenicsx4flow/pyvista4dolfinx/

(but it currently doesn’t support DG0 fields ?).

1 Like
5

Oh that’s pretty cool! So there are ways to visualize the data with matplotlib. That’s useful to know for someone like me that is very much still a beginner with all of this, but has several years of experience with matplotlib.

Yep, that result looks right. Thanks for your help!