I have a problem similar to Conditional statement in ufl form, where I am trying to define a conditional ufl function so that it can be used as a fem.Function. In this case, this is my analytic solution that I will compare the numerical solution to. I’m currently getting an error, and I would really appreciate some insight into what is going wrong, how I could make the code better, or if there is a different way of approaching this problem.
My conditional function is defined on the unit interval [0,1] and is
\begin{cases}
\alpha_1\mathcal{I}_1(x/\gamma) & 0 \leq x \leq r_1 \\
0.5 \mu_0 x\alpha_2 + \beta_2 / x & r_1\leq x\leq r_2 \\
\beta_3/x & r_2\leq x
\end{cases}
where \alpha_1, \alpha_2, \beta_2, \beta_3 are predetermined constants, and \mathcal{I}_1 is a modified bessel function of the first kind.
Please see below for a minimum working example and the error. Many thanks in advance!
Katie
P.s. Please note that this performed using the complex dolfinx kernel.
import numpy as np
from dolfinx import mesh, fem
import ufl
from mpi4py import MPI
domain = mesh.create_unit_interval(MPI.COMM_WORLD, 100)
V = fem.FunctionSpace(domain, ("Lagrange", 1))
x = ufl.SpatialCoordinate(domain)
analytic_A = fem.Function(V)
def analytic_rz(x):
# Simplified parameters for the problem
r1 = 0.1
r2 = 0.2
r3 = 0.3
mu0 = 4 * np.pi * 1e-7
C = 2e8
I = 100
gamma = np.sqrt(1j * C)
alpha2 = I # constants that are already known
beta1 = 0
alpha_ext = 0
alpha1 = (
mu0
* r2
* alpha2
/ (ufl.bessel_I(1, r1) + r1 * gamma * 0.5 * ( ufl.bessel_I(0, r1) + ufl.bessel_I(1, r1)))
)
beta2 = r1 * (alpha1 * ufl.bessel_I(1, r1) - 0.5 * mu0 * r1 * alpha2)
beta_ext = r2 * (beta2 / r2 - 0.5 * mu0 * r2 * alpha2)
A = ufl.conditional(ufl.And(0 <= x[0], x[0] <= r1), alpha1 * ufl.bessel_I(1, x[0] / gamma ), \
ufl.conditional(ufl.And(r1 <= x[0], x[0] <= r2), 0.5 * mu0 * x[0] * alpha2 + beta2 / x[0], beta_ext / x[0]))
return A
analytic_A.interpolate(analytic_rz)
TypeError Traceback (most recent call last)
File ~/anaconda3/envs/env3-10-complex/lib/python3.10/site-packages/dolfinx/fem/function.py:397, in Function.interpolate(self, u, cells, nmm_interpolation_data)
395 try:
396 # u is a Function or Expression (or pointer to one)
--> 397 _interpolate(u, cells)
398 except TypeError:
399 # u is callable
File ~/anaconda3/envs/env3-10-complex/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/env3-10-complex/lib/python3.10/site-packages/dolfinx/fem/function.py:373, in Function.interpolate.<locals>._interpolate(u, cells)
372 """Interpolate a cpp.fem.Function"""
--> 373 self._cpp_object.interpolate(u, cells, nmm_interpolation_data)
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
1. (self: dolfinx.cpp.fem.Function_complex128, f: numpy.ndarray[numpy.complex128], cells: numpy.ndarray[numpy.int32]) -> None
2. (self: dolfinx.cpp.fem.Function_complex128, u: dolfinx.cpp.fem.Function_complex128, cells: numpy.ndarray[numpy.int32], nmm_interpolation_data: Tuple[List[int], List[int], List[float], List[int]]) -> None
3. (self: dolfinx.cpp.fem.Function_complex128, expr: dolfinx::fem::Expression<std::complex<double>, double>, cells: numpy.ndarray[numpy.int32]) -> None
Invoked with: <dolfinx.cpp.fem.Function_complex128 object at 0x7f872bd772b0>, <function analytic_rz at 0x7f872bb07880>, 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), ((), (), (), ())
Did you forget to `#include <pybind11/stl.h>`? Or <pybind11/complex.h>,
<pybind11/functional.h>, <pybind11/chrono.h>, etc. Some automatic
conversions are optional and require extra headers to be included
when compiling your pybind11 module.
During handling of the above exception, another exception occurred:
ValueError Traceback (most recent call last)
Cell In[2], line 46
42 A = ufl.conditional(ufl.And(0 <= x[0], x[0] <= r1), alpha1 * ufl.bessel_I(1, x[0] / gamma ), \
43 ufl.conditional(ufl.And(r1 <= x[0], x[0] <= r2), 0.5 * mu0 * x[0] * alpha2 + beta2 / x[0], beta_ext / x[0]))
44 return A
---> 46 analytic_A.interpolate(analytic_rz)
File ~/anaconda3/envs/env3-10-complex/lib/python3.10/site-packages/dolfinx/fem/function.py:402, in Function.interpolate(self, u, cells, nmm_interpolation_data)
400 assert callable(u)
401 x = _cpp.fem.interpolation_coords(self._V.element, self._V.mesh.geometry, cells)
--> 402 self._cpp_object.interpolate(np.asarray(u(x), dtype=self.dtype), cells)
Cell In[2], line 42, in analytic_rz(x)
39 beta2 = r1 * (alpha1 * ufl.bessel_I(1, r1) - 0.5 * mu0 * r1 * alpha2)
41 beta_ext = r2 * (beta2 / r2 - 0.5 * mu0 * r2 * alpha2)
---> 42 A = ufl.conditional(ufl.And(0 <= x[0], x[0] <= r1), alpha1 * ufl.bessel_I(1, x[0] / gamma ), \
43 ufl.conditional(ufl.And(r1 <= x[0], x[0] <= r2), 0.5 * mu0 * x[0] * alpha2 + beta2 / x[0], beta_ext / x[0]))
44 return A
File ~/anaconda3/envs/env3-10-complex/lib/python3.10/site-packages/ufl/operators.py:500, in And(left, right)
498 def And(left, right):
499 """A boolean expression (left and right) for use with conditional."""
--> 500 return AndCondition(left, right)
File ~/anaconda3/envs/env3-10-complex/lib/python3.10/site-packages/ufl/conditional.py:209, in AndCondition.__init__(self, left, right)
207 def __init__(self, left, right):
208 """Initialise."""
--> 209 BinaryCondition.__init__(self, "&&", left, right)
File ~/anaconda3/envs/env3-10-complex/lib/python3.10/site-packages/ufl/conditional.py:50, in BinaryCondition.__init__(self, name, left, right)
48 def __init__(self, name, left, right):
49 """Initialise."""
---> 50 left = as_ufl(left)
51 right = as_ufl(right)
53 Condition.__init__(self, (left, right))
File ~/anaconda3/envs/env3-10-complex/lib/python3.10/site-packages/ufl/constantvalue.py:501, in as_ufl(expression)
499 return IntValue(expression)
500 else:
--> 501 raise ValueError(
502 f"Invalid type conversion: {expression} can not be converted to any UFL type.")
ValueError: Invalid type conversion: [ True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
True True True True True True True True] can not be converted to any UFL type.