How do I define a function with locally varying coefficients?

UFL
, ,
1

Hi everybody!
I am trying to define a function, depending on locally varyíng coefficients and on spatial coordinates, but so far without success. Here is my MWE based on a tutorial example “Defining subdomains for different materials”.

I am using FEniCSx v0.10.0.

from dolfinx import default_scalar_type
from dolfinx.fem import (
Constant,
Function,
functionspace,
)
from dolfinx.mesh import create_unit_square, locate_entities

from mpi4py import MPI

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
Q = functionspace(mesh, (“DG”, 0))

def Omega_0(x):
return x[1] <= 0.5

def Omega_1(x):
return x[1] >= 0.5

kappa = Function(Q)
cells_0 = locate_entities(mesh, mesh.topology.dim, Omega_0)
cells_1 = locate_entities(mesh, mesh.topology.dim, Omega_1)

myfunction = Function(Q)
myfunction.interpolate(lambda x: kappa*x[0])

But that does not work, the message is

UserWarning: Couldn’t map ‘f’ to a float, returning ufl object without evaluation.

I’ve figured out how to interpolate a function across subdomains and use that to create a global function. However, for my application, it would be great, if the coefficients depend on space and that the function can be interpolated with coefficients over the entire domain.

Does anyone have any idea how this might work? Is there a way? I’d really appreciate any feedback.

2

This is happening because you’re mixing numerical data (x[0] in your interpolation lambda) and symbolic algebra (kappa which is a finite element function defined on Q).

You can get around this by turning everything into a symbolic expression and interpolating. With some small modifications to your code:

from dolfinx import default_scalar_type
from dolfinx.fem import (
Constant,
Function,
functionspace,
Expression
)
from dolfinx.mesh import create_unit_square, locate_entities
import ufl

from mpi4py import MPI

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
Q = functionspace(mesh, ("DG", 0))

def Omega_0(x):
	return x[1] <= 0.5

def Omega_1(x):
	return x[1] >= 0.5

kappa = Function(Q)
cells_0 = locate_entities(mesh, mesh.topology.dim, Omega_0)
cells_1 = locate_entities(mesh, mesh.topology.dim, Omega_1)

myfunction = Function(Q)
x = ufl.SpatialCoordinate(mesh)
kappa.x.array[:] = 1.0
myfunction.interpolate(Expression(kappa*x[0], Q.element.interpolation_points))

where I set \kappa = 1 to just make something happen in the visualisation:

image
image821×670 13.8 KB

3

Hi nate!
Great, that works! Thank you very much for your help and your explanation.

Particularly because the code above wasn’t complete. Unfortunately, the assignment for kappa had been omitted:

kappa.x.array[cells_0] = np.full_like(cells_0, 1, dtype=default_scalar_type)
kappa.x.array[cells_1] = np.full_like(cells_1, 0.1, dtype=default_scalar_type)

Once again, thank you very much!