How to define initial conditions in mixed formulation using dolfinx?

dolfinx
1

Hello,

I would like to define initial conditions in preparation for a time dependent solve. In legacy fenics, this could be done via “y = interpolate(Constant((0.0, 1.0, 2.0)), ME)”. However, this does seem to work differently with the new dolfinx, since the following MWE results in a type error. Could you please advice what I have to change in order to make that work? Thank you.

from ufl import SpatialCoordinate, FiniteElement, MixedElement, TestFunctions, split
from dolfinx import mesh
from dolfinx.fem import FunctionSpace, Function, Constant
from mpi4py import MPI
import numpy

msh = mesh.create_interval(MPI.COMM_WORLD, 10, points=(0, 1))
x = SpatialCoordinate(msh)
CG1_elem = FiniteElement("CG", msh.ufl_cell(), 1)
ME_elem = MixedElement([CG1_elem, CG1_elem, CG1_elem])
ME = FunctionSpace(msh, ME_elem)
u_test, v_test, w_test = TestFunctions(ME)
y = Function(ME)
y.interpolate(Constant(msh, (0.0, 1.0, 2.0)), ME)
u, v, w = split(y)
2

See for example

3

Hello nate,
thank you for the examples. I am now able to use interpolation on a simple finite element function space. However, when I try to use the same concept on a mixed element function space, I get the following error: “Cannot get interpolation points - no Basix element available. Maybe this is a mixed element?” Here is the corresponding MWE.

from ufl import SpatialCoordinate, FiniteElement, MixedElement
from dolfinx import mesh
from dolfinx.fem import FunctionSpace, Function
from mpi4py import MPI
import numpy

def y_init(x):
    values = numpy.zeros((2, x.shape[1]))
    values[0] = 1.0
    values[1] = 2.0
    return values

msh = mesh.create_interval(MPI.COMM_WORLD, 10, points=(0, 1))

x = SpatialCoordinate(msh)

CG1_elem = FiniteElement("CG", msh.ufl_cell(), 1)
ME_elem = MixedElement([CG1_elem, CG1_elem])
ME = FunctionSpace(msh, ME_elem)

y = Function(ME)
y.interpolate(y_init)

Thank you for your help.

4

You should interpolate for each component of the problem:

def y0_init(x):
    values = numpy.zeros((1, x.shape[1]))
    values[0] = 1.0
    return values

def y1_init(x):
    values = numpy.zeros((1, x.shape[1]))
    values[0] = 2.0
    return values


msh = mesh.create_interval(MPI.COMM_WORLD, 10, points=(0, 1))

x = SpatialCoordinate(msh)

CG1_elem = FiniteElement("CG", msh.ufl_cell(), 1)
ME_elem = MixedElement([CG1_elem, CG1_elem])
ME = FunctionSpace(msh, ME_elem)

y = Function(ME)
y.sub(0).interpolate(y0_init)
y.sub(1).interpolate(y1_init)
1 Like