Using Project() with interpolated function

variational formulation
1

Hi everyone,

Currently withing my code, I have produced a interpolated function to represent a function changing with temperature for a 3D heat equation. I have produced this function using UnivariateSpline:

data = np.genfromtxt("text.txt")
k = data[:,0]
temp = data[:,1]
k_func = UnivariateSpline(temp, k, k = 2, ext = 3)

However I am now trying to project the value of this over my mesh using project():

k_space = FunctionSpace(mesh, 'CG',1)
k_new = project(k_f(T), kappa_space)

However I get the following error:

Traceback (most recent call last):
  File "test.py", line 213, in <module>
    k_new = project(k_f(T), kappa_space)
  File "/usr/lib/python3/dist-packages/scipy/interpolate/fitpack2.py", line 312, in __call__
    return fitpack.splev(x, self._eval_args, der=nu, ext=ext)
  File "/usr/lib/python3/dist-packages/scipy/interpolate/fitpack.py", line 368, in splev
    return _impl.splev(x, tck, der, ext)
  File "/usr/lib/python3/dist-packages/scipy/interpolate/_fitpack_impl.py", line 598, in splev
    y, ier = _fitpack._spl_(x, der, t, c, k, ext)
TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'

How can one get around this error whilst still being able to used this interpolated function?

A simplified version of my code can be found below:

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from ufl import as_tensor
from ufl import Index
import math
from scipy.optimize import curve_fit
from mpi4py import MPI
from scipy.interpolate import UnivariateSpline

mesh = UnitSquareMesh(10,10)

Space = FunctionSpace(mesh, 'P', 1)
T = project(Expression("x[0]+0.1*x[1]", degree=2), Space)

k_data = [30,40,50,60,70,80,90,100]
temp =[300,400,500,600,700,800,900,1000]

k_f = UnivariateSpline(temp, k_data, k = 1, ext = 3)

Space = FunctionSpace(mesh, 'P', 1)  
T = Function(Space)

kappa_space = FunctionSpace(mesh, 'CG',1)
k_new = project(k_f(T), kappa_space)

i = Index()
G = as_tensor(T.dx(i), (i))

q = -k_new*G
1 Like
2

Please create a minimal working code example that can be used to reproduce your error. Please follow the guidelines in Read before posting: How do I get my question answered?

3

This has now been added

4

Note that your function, k_f, takes an input that is an array of numbers, you are passing in a dolfin.Function.
What I suggest you do then, is the following:

k_new_array = k_f(T.vector().get_local())

Similarly, your project of k_f(T), an numpy array, is not an operation supported by dolfin. You need to either use setValuesLocal, as shown in
How to modify the value of a Function object? - #2 by dokken or

k_new = Function(Space)
k_new.vector()[:] = k_new_array
2 Likes
5

Thank you very much, this solves the problem.