How can I interpolate a UFL nabla_grad/grad to a function

dolfinx
1

Hello everyone and thank you in advance for your help.

I am trying to solve an equation and at the end, I have to use gradient (ufl.grad / ufl.nabla_grad) on the results by integrating the results.

I need to see the array after the “grad” or interpolate it on a function. However, I encounter this warning:

Couldn't map 'f' to a float, returning ufl object without evaluation.

here is a minimal example that produces the same warning. It is the same as demo_types.py and I only add ufl.grad(uh) to the end.

import numpy as np
import scipy.sparse
import scipy.sparse.linalg

import ufl
from dolfinx import fem, la, mesh, plot

from mpi4py import MPI

# -

# SciPy solvers do not support MPI, so all computations are performed on
# a single MPI rank

# +
comm = MPI.COMM_SELF
# -

# Create a mesh and function space

msh = mesh.create_rectangle(comm=comm, points=((0.0, 0.0), (2.0, 1.0)), n=(32, 16),
                            cell_type=mesh.CellType.triangle)
V = fem.FunctionSpace(msh, ("Lagrange", 1))

# Define a variational problem
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
fr = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
fc = ufl.sin(2 * np.pi * x[0]) + 10 * ufl.sin(4 * np.pi * x[1]) * 1j
gr = ufl.sin(5 * x[0])
gc = ufl.sin(5 * x[0]) * 1j
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = ufl.inner(fr + fc, v) * ufl.dx + ufl.inner(gr + gc, v) * ufl.ds


# In preparation for constructing Dirichlet boundary conditions, locate
# facets on the constrained boundary and the corresponding
# degrees-of-freedom
facets = mesh.locate_entities_boundary(msh, dim=1,
                                       marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),
                                                                      np.isclose(x[0], 2.0)))
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)


def solve(dtype=np.float32):
    """Solve the variational problem"""

    # Process forms. This will compile the forms for the requested type.
    a0 = fem.form(a, dtype=dtype)
    if np.issubdtype(dtype, np.complexfloating):
        L0 = fem.form(L, dtype=dtype)
    else:
        L0 = fem.form(ufl.replace(L, {fc: 0, gc: 0}), dtype=dtype)

    # Create a Dirichlet boundary condition
    bc = fem.dirichletbc(value=dtype(0), dofs=dofs, V=V)

    # Assemble forms
    A = fem.assemble_matrix(a0, [bc])
    A.finalize()

    b = fem.assemble_vector(L0)
    fem.apply_lifting(b.array, [a0], bcs=[[bc]])

    b.scatter_reverse(la.ScatterMode.add)
    fem.set_bc(b.array, [bc])

    # Create a Scipy sparse matrix that shares data with A
    As = scipy.sparse.csr_matrix((A.data, A.indices, A.indptr))

    # Solve the variational problem and return the solution
    uh = fem.Function(V, dtype=dtype)
    uh.x.array[:] = scipy.sparse.linalg.spsolve(As, b.array)
    um = fem.Function(V, dtype=dtype)
    um.interpolate(ufl.grad(uh))
    return um
2

See, e.g. dolfinx/demo_elasticity.py at main · FEniCS/dolfinx · GitHub

3

Thank you Nate. I could manage to solve the problem by the method introduced in the demo.

I’ve also noticed something else that maybe can help others (like me, new to dolfinx) when they want to solve PDE based on a mesh from gmsh.

when I created my mesh with gmsh in my python code from geo file and imported it by

msh, cell_markers, facet_markers = gmshio.model_to_mesh(model=model, comm=comm, rank=0)

I had a problem with the dimensions. The following commands resulted to:

print(msh.topology.dim, msh.geometry.dim). ==> (2,3)

and when I created VectorFunctionspace and a Function from it, the dimension of the Function became (3, )
The problem was solved by using “gdim=2” in the:

msh, cell_markers, facet_markers = gmshio.model_to_mesh(model=model, comm=comm, rank=0, gdim=2)

Then I could interpolate the result of ufl.nabla_grad(u) on my function without error.

I have also a question regarding Pyvista. How can I plot a function from “N1curl” space? Pyvista only supports CG and DG. In dolfin, I used to use dolfin.plot but I don’t know how to do it in dolfinx.

RuntimeError: Can only create meshes from continuous or discontinuous Lagrange spaces

Thanks

4

You need to interpolate a N1Curl function into a suitable vector DG space to visualize it with pyvista.
You can look at the compatibilities at the polynomial set definitions

2 Likes