Interpolating gradient of result in a particular direction and at a particular point

General
1

Hi everyone, I have 2 questions. Consider below the simple implementation of a Poisson equation.

from mpi4py import MPI
from dolfinx import mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

from dolfinx.fem import FunctionSpace
V = FunctionSpace(domain, ("Lagrange", 1))
from dolfinx import fem
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)

import numpy
# 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 = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)

import ufl
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
from dolfinx import default_scalar_type
f = fem.Constant(domain, default_scalar_type(-6))

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

Here are my 2 questions:

  • I have calculated the gradient of uh in X and Y directions by writing a code similar to what I saw here (Also attached the code below). I see that this is particularly useful when making calculations using the gradients. But, is there an easier way of extracting the gradient in X & Y direction? My objective is to be able to map and find the gradient on each direction to a point in the mesh.
  • How to interpolate the X and Y gradients at a random point (let’s say (0.5, 0.5)) on the geometry? I am looking for a way similar to the function ‘evaluateGradient’ in MATLAB.

MWE for gradient in X-direction:

x = SpatialCoordinate(domain)
x_grad = inner(as_vector((1, 0)), grad(uh))
W = fem.FunctionSpace(domain, ("DQ", 1))
expr = fem.Expression(x_grad, W.element.interpolation_points())
w = fem.Function(W)
w.interpolate(expr)

Thank you!

2

This can be simplified to
x_grad = grad(uh)[0]
or simply
x_grad = uh.dx(0)
Im not sure what

Is supposed to do in the code above, as it is not used.

For the last part of your question, see:
https://jsdokken.com/dolfinx-tutorial/chapter1/membrane_code.html#making-curve-plots-throughout-the-domain
which explains how to evaluate a function at any set of points, thus after interpolating the gradient into an appropriate space, you can use the eval-functionality.

There are other ways of Getting the gradient at a given point in a slightly more efficient fashion, but as Im currently not at my computer, I dont have time to sketch out another version.

3

Thank you @dokken . This helps!

Yes, the x = SpatialCoordinate(domain) line is not needed. I removed it!

If it is possible, can you share the more efficient way of doing this? I am asking because I utilize this evaluation in an implementation which solves my model in an iterative fashion. When doing it in MATLAB, each time step takes 3-4 minutes to come up with the solution and I do it for at least 100 time steps. Any slight improvement in efficiency would make a difference!

Thank you!

4

Dear @appana_3,
Could you create a minimal script with timings that illustrate what part of your code is time-consuming?
This would help me help you address any performance issues.