Map solution to another function space

dolfinx
1

I used to be able to map a solution function from one function space to another (for exporting to paraview) using function

def interpolate(V, f):
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
    a_p = ufl.inner(u, v) * ufl.dx
    L_p = ufl.inner(f, v) * ufl.dx
    return fem.petsc.LinearProblem(a_p, L_p).solve()

V = fem.functionspace(msh, mixed_element([N1curl, H1]))
Vproj2d = fem.functionspace(msh, ("Lagrange", 2))
projection = interpolate(Vproj2d, solution)

In newer DolfinX version, this fails with error

ValueError: Shapes do not match: <Coefficient id=140714964273552> and <Argument id=140715103632448>
Backtrace 
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[11], line 13
     10     continue
     11 print(f"eigenvalue: {eigenval:.12f};\t"
     12       f"int(psi): {dolfinx.fem.assemble_scalar(fem.form(E2 * ufl.dx)):.12e}")
---> 13 func = interpolate(Vproj2d, E1)
     14 func.name = f"E1-{j:03d}-{eigenval:.4f}"
     15 xdmf.write_function(func)

Cell In[2], line 4, in interpolate(V, f)
      2 u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
      3 a_p = ufl.inner(u, v) * ufl.dx
----> 4 L_p = ufl.inner(f, v) * ufl.dx
      5 return fem.petsc.LinearProblem(a_p, L_p).solve()

File /opt/conda/lib/python3.11/site-packages/ufl/operators.py:213, in inner(a, b)
    211 if a.ufl_shape == () and b.ufl_shape == ():
    212     return a * Conj(b)
--> 213 return Inner(a, b)

File /opt/conda/lib/python3.11/site-packages/ufl/tensoralgebra.py:166, in Inner.__new__(cls, a, b)
    164 ash, bsh = a.ufl_shape, b.ufl_shape
    165 if ash != bsh:
--> 166     raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
    168 # Simplification
    169 if isinstance(a, Zero) or isinstance(b, Zero):

ValueError: Shapes do not match: <Coefficient id=140714964273552> and <Argument id=140715103632448>

How to do this in newer versions of DolfinX?

Full code 
import dolfinx
from dolfinx import io, fem, mesh
from basix.ufl import element, mixed_element
import ufl
import numpy as np
from mpi4py import MPI
from slepc4py import SLEPc
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import assemble_matrix as petsc_assemble_matrix

def interpolate(V, f):
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
    a_p = ufl.inner(u, v) * ufl.dx
    L_p = ufl.inner(f, v) * ufl.dx
    return fem.petsc.LinearProblem(a_p, L_p).solve()

msh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD,
    [np.array([0, 0]), np.array([np.pi, np.pi])], [40, 40],
    cell_type=dolfinx.cpp.mesh.CellType.triangle)
mdim = msh.topology.dim - 1
msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)

N1curl = element("Nedelec 1st kind H(curl)", msh.basix_cell(), 1)
H1 = element("Lagrange", msh.basix_cell(), 1)
V = fem.functionspace(msh, mixed_element([N1curl, H1]))
Vproj2d = fem.functionspace(msh, ("Lagrange", 2))

e, p = ufl.TrialFunctions(V)
psi, q = ufl.TestFunctions(V)

a = ufl.inner(ufl.curl(e), ufl.curl(psi)) * ufl.dx + ufl.inner(ufl.grad(p), psi) * ufl.dx
a += ufl.inner(e, ufl.grad(q)) * ufl.dx
b = ufl.inner(e, psi) * ufl.dx

bc_facets = dolfinx.mesh.exterior_facet_indices(msh.topology)
bc_dofs = fem.locate_dofs_topological(V, mdim, bc_facets)

u_bc = fem.Function(V)
u_bc.x.array[:] = ScalarType(0)
bc = fem.dirichletbc(u_bc, bc_dofs)

A = petsc_assemble_matrix(fem.form(a), bcs=[bc], diagonal=0.0)
A.assemble()
B = petsc_assemble_matrix(fem.form(b), bcs=[bc])
B.assemble()

eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, B)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
st = eps.getST()
st.setType(SLEPc.ST.Type.SINVERT)
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
eps.setTarget(10)
eps.setDimensions(20)
eps.solve()

vals = [(i, eps.getEigenvalue(i)) for i in range(eps.getConverged())]
vals.sort(key=lambda x: x[1].real)

j = 0
E = fem.Function(V)
with io.XDMFFile(MPI.COMM_WORLD, "out.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh)
    for i, _ in vals:
        E1, E2 = E.split()
        eigenval = eps.getEigenpair(i, E.vector).real
        penalty = np.abs(fem.assemble_scalar(fem.form(E2 * ufl.dx)))
        if not np.isclose(penalty, 0):
            continue
        print(f"eigenvalue: {eigenval:.12f};\t"
              f"int(psi): {dolfinx.fem.assemble_scalar(fem.form(E2 * ufl.dx)):.12e}")
        func = interpolate(Vproj2d, E1)
        func.name = f"E1-{j:03d}-{eigenval:.4f}"
        xdmf.write_function(func)
        func = interpolate(Vproj2d, E1)
        func.name = f"E2-{j:03d}-{eigenval:.4f}"
        xdmf.write_function(func)
        j += 1
2

The first function space in your mixed element is a vector space:

Which is not compatible with Vproj2d as it is a scalar space.
You can fix this by adding:

Vproj2d = fem.functionspace(msh, ("Lagrange", 2, (msh.geometry.dim, )))

However, note that a second order continuous Lagrange element is not suitable to represent Nedelec spaces, and one should use an appropriate DG space and VTKFile or VTXWriter for proper visualization.

3

Thanks! I meant to use Lagrange in 2D space, not 2nd order.

Now, paraview throws segfault. It’s always something…