Interpolating step functions in fenicsx

General
1

Hi all,

I’m trying to interpolate a step function as part of some fenicsx code. When I do, I get an error, so I’m clearly approaching this the wrong way.

I’m giving a MWE below: can you tell me what’s wrong?

Thanks in advance!

from mpi4py import MPI
import dolfinx as d
import ufl
import matplotlib.pyplot as plt
import numpy as np

msh = d.mesh.create_interval( comm=MPI.COMM_WORLD,
                                    nx=100, points=(0, 10) )
Pn = ufl.FiniteElement( 'CG', msh.ufl_cell(), 4 )
S = d.fem.FunctionSpace( msh, Pn )
f = d.fem.Function(S)

def profile( r, rho0, r_step ):
    if r < r_step:
        return rho0
    else:
        return 0.

# some calculations to determine rho0 and r_step...
# here I assign some random values, just to make a MWE

prof = lambda r : profile( r[0], 10, 1 )

f.interpolate(prof)

The error I get is:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
File ~/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/function.py:370, in Function.interpolate(self, u, cells)
    368 try:
    369     # u is a Function or Expression (or pointer to one)
--> 370     _interpolate(u, cells)
    371 except TypeError:
    372     # u is callable

File ~/anaconda3/envs/fenicsx/lib/python3.10/functools.py:889, in singledispatch.<locals>.wrapper(*args, **kw)
    886     raise TypeError(f'{funcname} requires at least '
    887                     '1 positional argument')
--> 889 return dispatch(args[0].__class__)(*args, **kw)

File ~/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/function.py:346, in Function.interpolate.<locals>._interpolate(u, cells)
    345 """Interpolate a cpp.fem.Function"""
--> 346 self._cpp_object.interpolate(u, cells)

TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfinx.cpp.fem.Function_float64, f: numpy.ndarray[numpy.float64], cells: numpy.ndarray[numpy.int32]) -> None
    2. (self: dolfinx.cpp.fem.Function_float64, u: dolfinx.cpp.fem.Function_float64, cells: numpy.ndarray[numpy.int32]) -> None
    3. (self: dolfinx.cpp.fem.Function_float64, expr: dolfinx::fem::Expression<double>, cells: numpy.ndarray[numpy.int32]) -> None

Invoked with: <dolfinx.cpp.fem.Function_float64 object at 0x1665b4bb0>, <function <lambda> at 0x15adb0790>, array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
       51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,
       68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84,
       85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99],
      dtype=int32)

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
/path_to_file/test.ipynb Cell 13 line 2
     23 # some calculations to determine rho0 and r_step...
     24 # here I assign some random values, just to make a MWE
     26 prof = lambda r : profile( r[0], 10, 1 )
---> 28 f.interpolate(prof)

File ~/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/function.py:377, in Function.interpolate(self, u, cells)
    373 assert callable(u)
    374 x = _cpp.fem.interpolation_coords(
    375     self._V.element, self._V.mesh, cells)
    376 self._cpp_object.interpolate(
--> 377     np.asarray(u(x), dtype=self.dtype), cells)

/path_to_file/test.ipynb Cell 13 line 2
     21         return 0.
     23 # some calculations to determine rho0 and r_step...
     24 # here I assign some random values, just to make a MWE
---> 26 prof = lambda r : profile( r[0], 10, 1 )
     28 f.interpolate(prof)

/path_to_file/test.ipynb Cell 13 line 1
     17 def profile( r, rho0, r_step ):
---> 18     if r < r_step:
     19         return rho0
     20     else:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
2

The interpolation function in DOLFINx is vectorized. It means that it takes in multiple points at once, i.e. the input r is a (num_points, 3) dimensional array. Thus, you would either have to vectorize your code, i.e.

def profile( r, rho0, r_step ):
    return rho0* (r < r_step)

which for an simple a:

 a = np.arange(10)
profile(a, 0.1, 3)

returns

Out[10]: array([0.1, 0.1, 0.1, 0. , 0. , 0. , 0. , 0. , 0. , 0. ])
1 Like
3

Thanks! I had sensed the function was being vectorised, and I had tried, naïvely, with f.interpolate( np.vectorize(prof) ), but it hadn’t worked. That makes sense. Thank you!