How to return "ufl.algebra" instead of "numpy.ndarray"?

dolfinx
1

Greetings,

I need some help with “ufl.algebra” class. The following code introduces two functions:

import dolfinx
import numpy as np
from mpi4py import MPI
import scipy.interpolate as sp

mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [1, 1])
V = dolfinx.FunctionSpace(mesh, ("CG", 1))
uh = dolfinx.Function(V)
# Fill with some values
uh.interpolate(lambda x: x[0]+1)

def polynom(var):
    return 0.00004*var + 0.02 

b = polynom(uh)
print(type(b))
print(b)

def interpolation(var):
    uh_extracted_values = np.array(var.vector.array) 
    data_points = sp.interp1d([0,1000],[0.02,0.06])
    return data_points(uh_extracted_values)

b = interpolation(uh)
print(type(b))
print(b)

output:

<class 'ufl.algebra.Sum'>
0.02 + 4e-05 * f_2
<class 'numpy.ndarray'>
[0.02004 0.02004 0.02012 0.02012]

Does anyone know how to return “ufl.algebra” class for the “interpolation” function? I cannot use the class “numpy.ndarray” further in my simulations since I encounter the following error:

Traceback (most recent call last):
  File "/root/weldingfenicsx/test_nonlinearity.py", line 67, in <module>
    problem = dolfinx.fem.NonlinearProblem(F, uh)
  File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/problem.py", line 171, in __init__
    self._L = fem.form.Form(F, form_compiler_parameters=form_compiler_parameters, jit_parameters=jit_parameters)
  File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/form.py", line 34, in __init__
    sd = form.subdomain_data()
AttributeError: 'numpy.ndarray' object has no attribute 'subdomain_data'

It seems like passing dolfinx.Function as an input for a function makes an “ufl” output automatically, as shown below:

import dolfinx
import numpy as np
from mpi4py import MPI
import scipy.interpolate as sp

mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [1, 1])
V = dolfinx.FunctionSpace(mesh, ("CG", 1))
uh = dolfinx.Function(V)
# Fill with some values
uh.interpolate(lambda x: x[0]+1)

some_numbers = np.array([1, 1, 3, 3])

def polynom(var):
    return 0.00004*var + 0.02 

#dolfinx function
b = polynom(uh)
print(type(b))
print(b)

#numpy array
b = polynom(some_numbers)
print(type(b))
print(b)

output:

<class 'ufl.algebra.Sum'>
0.02 + 4e-05 * f_2
<class 'numpy.ndarray'>
[0.02004 0.02004 0.02012 0.02012]

The problem is that I cannot pass dolfinx.Function into interp1d:

import dolfinx
import numpy as np
from mpi4py import MPI
import scipy.interpolate as sp

mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [1, 1])
V = dolfinx.FunctionSpace(mesh, ("CG", 1))
uh = dolfinx.Function(V)
# Fill with some values
uh.interpolate(lambda x: x[0]+1)

def interpolation(var):
    data_points = sp.interp1d([0,1000],[0.02,0.06])
    return data_points(var)

b = interpolation(uh)
print(type(b))
print(b)

output:

Traceback (most recent call last):
  File "/root/weldingfenicsx/for_forum.py", line 16, in <module>
    b = interpolation(uh)
  File "/root/weldingfenicsx/for_forum.py", line 14, in interpolation
    return data_points(var)
  File "/usr/lib/python3/dist-packages/scipy/interpolate/polyint.py", line 78, in __call__
    x, x_shape = self._prepare_x(x)
  File "/usr/lib/python3/dist-packages/scipy/interpolate/polyint.py", line 90, in _prepare_x
    x = _asarray_validated(x, check_finite=False, as_inexact=True)
  File "/usr/lib/python3/dist-packages/scipy/_lib/_util.py", line 242, in _asarray_validated
    raise ValueError('object arrays are not supported')
ValueError: object arrays are not supported

Can anyone help me out?

Thanks

Sincerely,
Anton

2

I am not sure what you want to do with this code after the b=interpolation(uh), but you can avoid the error message with the following code:

import dolfinx
import numpy as np
from mpi4py import MPI
import scipy.interpolate as sp

mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [1, 1])
V = dolfinx.FunctionSpace(mesh, ("CG", 1))
uh = dolfinx.Function(V)
# Fill with some values
uh.interpolate(lambda x: x[0]+1)

def interpolation(var):
    data_points = sp.interp1d([0,1000],[0.02,0.06])
    return data_points(var.vector.array)

b = interpolation(uh)
print(type(b))
print(b)
2 Likes
3

Hi Jorgen,

in general I am trying to define material properties (as a function of solution) using data points instead of polynomials in dolfinx, see here. The solver works fine with the “polynom” function (which has an output type: ufl.algebra). However if I use the “interpolation” function to define properties (output type: numpy.ndarray) the solver gives the following error: “AttributeError: ‘numpy.ndarray’ object has no attribute ‘subdomain_data’”. Therefore I though that if the “interpolation” function returns exactly the same “ufl.algebra” data type it could work.

MWE:

import dolfinx
import numpy as np
import ufl
from mpi4py import MPI
import scipy.interpolate as sp

mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [1, 1])
V = dolfinx.FunctionSpace(mesh, ("CG", 1))
uh = dolfinx.Function(V)
u_n = dolfinx.Function(V)

# def polynom(u):
#     return 0.00004*u + 0.02

def interpolation(var):
    data_points = sp.interp1d([0,1000],[0.02,0.06])
    return data_points(var.vector.array)

v = ufl.TestFunction(V)
f = dolfinx.Constant(mesh, 0)

# F = polynom(uh)*ufl.dot(ufl.grad(uh), ufl.grad(v))*ufl.dx
F = interpolation(uh)*ufl.dot(ufl.grad(uh), ufl.grad(v))*ufl.dx

problem = dolfinx.fem.NonlinearProblem(F, uh)
solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
n, converged = solver.solve(uh)

Can I somehow adapt the “interpolation” function, such that it could be used for the existing nonlinear solver?

I could add:

q_ = dolfinx.Function(V)
q_.vector.array[:] = interpolation(uh)
F = q_*ufl.dot(ufl.grad(uh), ufl.grad(v))*ufl.dx

which doesn’t arise any error, however does not update material properties after each iteration as I showed here