Dolfinx discontinous expression

dolfinx
1

Hi all, I want to define a discontinous expression on two subdomain. The expression equals to 0 and 1 in the two subdomain respectively. How to realize it using dolfinx? Thank you in advance.

2

Hi, is this a expression that can be expressed geometrically, like Domain a 0.5<=x[1], Domain 2 x[1]=<0.5, or is it defined through a mesh function?

3

You could for instance do the following for geometrical search:

import dolfinx
import dolfinx.io

mesh = dolfinx.UnitSquareMesh(dolfinx.MPI.comm_world, 9, 9)
V = dolfinx.FunctionSpace(mesh, ("DG", 0))
v = dolfinx.Function(V)
x = V.tabulate_dof_coordinates()
for i in range(x.shape[0]):
    midpoint = x[i,:]
    if midpoint[0]> 0.5 and midpoint[1]>0.25:
        v.vector.setValueLocal(i, 2)
    else:
        v.vector.setValueLocal(i, 1)

dolfinx.io.XDMFFile(dolfinx.MPI.comm_world, "DG.xdmf").write_checkpoint(v, "v", 0)

and this if the data comes from a MeshFunction

v2 = dolfinx.Function(V)
mf = dolfinx.MeshFunction("size_t", mesh, mesh.topology.dim, 0)
mf.mark(lambda x: numpy.logical_and(x[0] < 0.5, x[1]>0.2), 2)
for i in range(x.shape[0]):
    if mf.values[i] == 2:
        v.vector.setValueLocal(i, 2)
    else:
        v.vector.setValueLocal(i, 1)

dolfinx.io.XDMFFile(dolfinx.MPI.comm_world, "DG_mf.xdmf").write_checkpoint(v, "v", 0)
3 Likes
4

dokken, it’s the answer I want. Many thanks.

5

Dolfinx is more pythonic and convenient than dolfin. I like it.

6

That is one of the goals of dolfinx. For instance, you can use numba to create custom assemblers, use normal Python functions (lambda functions) to locate dofs, integrate expressions etc.

7

Numba could be used in this case to accelerate the loops.

2 Likes
8

Hi everyone, sorry if I re-open this but I was wondering how to avoid the for loop indicated by @dokken (just to save some lines).

I would like something like:

foo = dolfinx.fem.Function(V)
foo.interpolate(lambda x: 1. if x[0] < 0.5 else 0)

But this does not work. The only way I found to avoid the for loop is to do the following:

x = ufl.SpatialCoordinate(mesh)
foo_ufl_exp = ufl.conditional(ufl.lt(x[0], 0.5), 1., 0)  # stands for 'x[0] < 0.5 ? 1 : 0'
foo_exp = dolfinx.fem.Expression(foo_ufl_exp,
                                 function_space.element.interpolation_points)
# interpolate
c_0.interpolate(c_0_exp)

But this look even less clear than the for loop. Anything I am missing?

9

See numpy.where.

1 Like
10

foo.interpolate(lambda x: np.where(x[0] < 0.5, 1, 0)) does the job indeed.

Thank you very much @nate