# Problem interpolating mixed function dolfinx

Hi,

I’m facing a problem imposing initial conditions in a mixed function space. When I open the files in ParaView, the values of some components were mixed in some nodes. I noticed that, without the VectorElement with dimension 3 and the TensorElement, the interpolation is ok.

``````from ufl import VectorElement, FiniteElement, TensorElement, MixedElement
import dolfinx
from dolfinx import RectangleMesh, FunctionSpace, Function
from dolfinx.io import XDMFFile
from mpi4py import MPI
import numpy as np

# Mesh
mesh = RectangleMesh(
MPI.COMM_WORLD,
[np.array([0, 0, 0]), np.array([1, 1, 0])], [32, 32],
dolfinx.cpp.mesh.CellType.triangle, dolfinx.cpp.mesh.GhostMode.none)

# Define Function Space
V = VectorElement("Lagrange", mesh.ufl_cell(), 1)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
S = VectorElement("Lagrange", mesh.ufl_cell(), 1, 3)
L = TensorElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([V, Q, S, L]))

# Define function
w = Function(W)

# Initial conditions
def w_init(x):
values = np.zeros((10, x.shape[1]))
values[0] = 1.0
values[1] = 2.0
values[2] = 3.0
values[3] = 4.0
values[4] = 5.0
values[5] = 6.0
values[6] = 7.0
values[7] = 8.0
values[8] = 9.0
values[9] = 10.0
return values

w.interpolate(w_init)

# Files
u_file = XDMFFile(MPI.COMM_WORLD, "test/velocity.xdmf", "w")
u_file.write_mesh(mesh)
u_file.write_function(w.split()[0])

p_file = XDMFFile(MPI.COMM_WORLD, "test/pressure.xdmf", "w")
p_file.write_mesh(mesh)
p_file.write_function(w.split()[1])

A11_file = XDMFFile(MPI.COMM_WORLD, "test/A11.xdmf", "w")
A11_file.write_mesh(mesh)
A11_file.write_function(w.split()[2].split()[0])

A12_file = XDMFFile(MPI.COMM_WORLD, "test/A12.xdmf", "w")
A12_file.write_mesh(mesh)
A12_file.write_function(w.split()[2].split()[1])

A22_file = XDMFFile(MPI.COMM_WORLD, "test/A22.xdmf", "w")
A22_file.write_mesh(mesh)
A22_file.write_function(w.split()[2].split()[2])

G11_file = XDMFFile(MPI.COMM_WORLD, "test/G11.xdmf", "w")
G11_file.write_mesh(mesh)
G11_file.write_function(w.split()[3].split()[0])

G12_file = XDMFFile(MPI.COMM_WORLD, "test/G12.xdmf", "w")
G12_file.write_mesh(mesh)
G12_file.write_function(w.split()[3].split()[1])

G21_file = XDMFFile(MPI.COMM_WORLD, "test/G21.xdmf", "w")
G21_file.write_mesh(mesh)
G21_file.write_function(w.split()[3].split()[2])

G22_file = XDMFFile(MPI.COMM_WORLD, "test/G22.xdmf", "w")
G22_file.write_mesh(mesh)
G22_file.write_function(w.split()[3].split()[3])

u_file.close()
p_file.close()
A11_file.close()
A12_file.close()
A22_file.close()
G11_file.close()
G12_file.close()
G21_file.close()
G22_file.close()
``````

These are the components x and y of the file velocity.xdmf, which are supposed to be equal to 1.0 and 2.0, respectively.

1 Like

Hello,

This is due to a bug in dofmap creation when using a `MixedElement`.

I’ve pushed a fix in this branch: https://github.com/FEniCS/dolfinx/tree/mscroggs/mixed_bugfix. I’ll open a pull request once I’ve added a unit test to make sure we don’t break this again.

Edit: Sorry, spoke too soon; it might not be fixed yet as tests are failing it that branch. I’ll update here once I make progress

1 Like

Thank you, for your reply! Is there a workaround to impose these initial conditions? I tried to use w.sub(0).vector.set(1.0) for example, but it change all de function values.

I have now fixed my fix: https://github.com/FEniCS/dolfinx/pull/1137. Hopefully this will be merged soon

I don’t think there’s a workaround in the current version, as any MixedElement containing anything other than a FiniteElement will assign the dof numbers incorrectly.

3 Likes

Hi,

I’m trying to impose a global constraint in my domain via Lagrange multiplier. For that, I have created a mixed function space composed of “CG” and “Real” elements. When I try to impose a initial condition interpolating the function created using the mixed space, I receive the following message:

``````Loguru caught a signal: SIGSEGV
Segmentation fault (core dumped)
``````

If I change “Real” for “DG” the code runs normally and the interpolation is ok. I decided to revive this topic because I think it could be some problem with the mixed function space. Here is a minimal working example:

``````import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
from dolfinx.cpp.mesh import CellType
from dolfinx import (RectangleMesh, Function, FunctionSpace)
from ufl import (FiniteElement, MixedElement, split)

mesh = RectangleMesh(
MPI.COMM_WORLD,
[np.array([0, -0.5, 0]), np.array([10, 0.5, 0])], [5, 5],

N = FiniteElement("CG", mesh.ufl_cell(), 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([N, R]))

w = Function(W)

def w_init(x):
values = np.zeros((2, x.shape[1]))
values[0] = 0.0
values[1] = 0.0
return values

w.interpolate(w_init)``````

I have managed to impose the initial conditions using projection dolfiny, the result is not as smooth as interpolating but it works.

``````def project(v, target_func, bcs=[]):
# Ensure we have a mesh and attach to measure
V = target_func.function_space
dx = ufl.dx(V.mesh)

# Define variational problem for projection
w = ufl.TestFunction(V)
Pv = ufl.TrialFunction(V)
a = ufl.inner(Pv, w) * dx
L = ufl.inner(v, w) * dx

# Assemble linear system
A = dolfinx.fem.assemble_matrix(a, bcs)
A.assemble()
b = dolfinx.fem.assemble_vector(L)
dolfinx.fem.apply_lifting(b, [a], [bcs])
dolfinx.fem.set_bc(b, bcs)

solve(A, target_func.vector, b)

w_init = dolfinx.Constant(mesh, (1.0, 2.0))
project(w_init, w)``````

Thanks for this minimal example.

It look like the Real space is broken in dolfinx: it should have one DOF in total, but appears to have one DOF per cell. I’ve opened an issue in the dolfinx repo and will try to get to the bottom of this.

1 Like

Hi, @mscroggs. I took some time trying to figure out a workaround for interpolation. After “solving” this I noticed that the Real space is not working properly and start a new topic link. Thank you for take part of your time to work on this issue!

I’ve posted a link to the issue on that thread to so anyone who finds that thread and not this one also sees it.

1 Like

Dear @mscroggs ,

I’m currently facing the same problem using
fenics-basix 0.7.0
fenics-dolfinx 0.7.3
fenics-ffcx 0.7.0
fenics-ufl 2023.2.0

from anaconda on Ubuntu 22.04.3 LTS, 22.04 jammy.

Since the OP uses an older fenicsx version, I created my own MRE. However, the interpolation on a mixed function space still seems to produce wrong results.

In particular, I’m using a mixed space with two vector spaces and a scalar space. I’m failing to interpolate a vector function `ic(x)` to the third subspace as well as a constant vector field to the first subspace. Furthermore, interpolating function `ic(x)` to another function `ic_values` yields bugs.

Here’s my example:

``````from dolfinx import io
from dolfinx import fem
from dolfinx.mesh import create_rectangle
from mpi4py import MPI
import numpy as np
import ufl

mesh = create_rectangle(
MPI.COMM_WORLD,
[np.array([-20, -20]), np.array([20, 20])], [32, 32])

xdmfout = io.XDMFFile(MPI.COMM_WORLD, 'output.xdmf', 'w')
xdmfout.write_mesh(mesh)
xdmfout.close()

coor = mesh.geometry.x
x0min=np.min(coor[:,0])

Velem = ufl.VectorElement('Lagrange', mesh.ufl_cell(), degree=1, dim=2)
Felem = ufl.FiniteElement('Lagrange', mesh.ufl_cell(), 1)

V = fem.functionspace(mesh, Velem )
F = fem.functionspace(mesh, Felem )
W = fem.FunctionSpace(mesh, ufl.MixedElement([Velem,Felem,Velem]))

w = fem.Function(W)

def ic(x):
'''vector that points downwards in 1/4 of the domain, upwards elsewhere'''
values = np.zeros((2, x.shape[1]))
values[0,:] = 0.0
values[1,:] = np.sign(x[0]-x0min*3/4)
return values

ic_values = fem.Function(V, name='ic_values')
ic_values.interpolate(ic)

w.sub(0).interpolate(fem.Expression(fem.Constant(mesh, [1.0, 2.0]), W.sub(0).element.interpolation_points()))
w.sub(1).interpolate(fem.Expression(fem.Constant(mesh, 3.0), W.sub(1).element.interpolation_points())) # this one works
w.sub(2).interpolate(ic)

u, phi, p = w.split()
u.name = 'u'
phi.name = 'phi'
p.name = 'p'

xdmfout = io.XDMFFile(MPI.COMM_WORLD, 'output.xdmf', 'a')

xdmfout.write_function(u, 0)
xdmfout.write_function(phi, 0)
xdmfout.write_function(p, 0)
xdmfout.write_function(ic_values, 0)

xdmfout.close()
``````

Results:
The interpolation of a constant vector field to `u`, the first subspace, shows similar bugs as in the OP.
The interpolation of a scalar field to the second subspace looks good.
The interpolation of a vector function to `p`, the third subspace, shows that the values were assigned falsely, namely the x and y component switched.
Interpolation to the independent function `ic_values` looks good, except the undesired corner at the top boundary.

I’d appreciate any help on this topic.
Best regards.

You should collapse the subfunctions prior to writing to file, i.e.

``````from dolfinx import io
from dolfinx import fem
from dolfinx.mesh import create_rectangle
from mpi4py import MPI
import numpy as np
import ufl

mesh = create_rectangle(
MPI.COMM_WORLD,
[np.array([-20, -20]), np.array([20, 20])], [32, 32])

coor = mesh.geometry.x
x0min=np.min(coor[:,0])

Velem = ufl.VectorElement('Lagrange', mesh.ufl_cell(), degree=1, dim=2)
Felem = ufl.FiniteElement('Lagrange', mesh.ufl_cell(), 1)

V = fem.functionspace(mesh, Velem )
F = fem.functionspace(mesh, Felem )
W = fem.FunctionSpace(mesh, ufl.MixedElement([Velem,Felem,Velem]))

w = fem.Function(W)

def ic(x):
'''vector that points downwards in 1/4 of the domain, upwards elsewhere'''
values = np.zeros((2, x.shape[1]))
values[0,:] = 0.0
values[1,:] = np.sign(x[0]-x0min*3/4)
return values

ic_values = fem.Function(V, name='ic_values')
ic_values.interpolate(ic)
w.sub(0).interpolate(fem.Expression(fem.Constant(mesh, [1.0, 2.0]), W.sub(0).element.interpolation_points()))
w.sub(1).interpolate(fem.Expression(fem.Constant(mesh, 3.0), W.sub(1).element.interpolation_points())) # this one works
w.sub(2).interpolate(ic)

u, phi, p = w.split()
u_c = u.collapse()
u_c.name = "u"
phi_c = phi.collapse()
phi_c.name = "phi"
p_c = p.collapse()
p_c.name = "p"
xdmfout = io.XDMFFile(MPI.COMM_WORLD, 'output.xdmf', 'a')
xdmfout.write_mesh(mesh)
xdmfout.write_function(u_c, 0)
xdmfout.write_function(phi_c, 0)
xdmfout.write_function(p_c, 0)
xdmfout.write_function(ic_values, 0)

xdmfout.close()
``````

Great, that already solved the problem. Thank you very much.