Problem interpolating mixed function dolfinx

fforteneto · September 23, 2020, 1:17am

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.

mscroggs · September 23, 2020, 10:41am

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

fforteneto · September 23, 2020, 11:48am

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.

mscroggs · September 23, 2020, 3:24pm

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.

fforteneto · October 27, 2020, 6:22pm

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],
    CellType.quadrilateral, dolfinx.cpp.mesh.GhostMode.none)

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)

fforteneto · November 11, 2020, 6:10pm

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])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    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)

mscroggs · November 12, 2020, 9:28am

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.

fforteneto · November 12, 2020, 11:32am

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!

mscroggs · November 12, 2020, 11:54am

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.

mcjbo · April 15, 2024, 1:47pm

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.

dokken · April 15, 2024, 2:30pm

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()

mcjbo · April 15, 2024, 3:14pm

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

Topic		Replies	Views
Interpolate Constant as Initial Condition in mixed element space Errors	6	977	December 12, 2022
Set initial conditions for tensor-element dolfinx	6	781	July 17, 2022
Interpolate into complex-valued vector function dolfinx	3	784	September 2, 2022
Initial Conditions and Function update on time in mixed Element dolfinx	2	275	December 5, 2023
Quadrature function space dolfinx	8	1519	December 19, 2022

Problem interpolating mixed function dolfinx

Related Topics