Using scipy special functions

1

Hello everyone,

I’m trying to use some special functions from scipy, namely, the first and second kind elliptic integrals which, as far as I know, are not implemented in ufl.

Specifically I want to project an analytical function that depends in this functions into my FunctionSpace.

Thanks in advance.

2

First of all, you need to specify if you use legacy DOLFIN or DOLFINx.

Secondly, a minimal (not working) example of what you want to achieve, on a standard geometry (cube or square) would make it easier for people to help you

3

True, sorry for the inacuracy, I’m using legacy DOLFIN.

Take this code snippet as an example:

import dolfin as df
import numpy as np
import ufl
from scipy.special import ellipk, ellipe
import matplotlib.pyplot as plt

# Define the mesh and function spaces

mesh = df.UnitSquareMesh(50, 50)
cg_degree = 2
V = df.FunctionSpace(mesh, 'CG', cg_degree)

def f1(z=0., x=0., lib=np):
      '''
      Returns function f at a given (set of) point(s).
      '''
      f1= z*0. # allocate and zero
      f1 = x**2 + z**2
      return f1

def f2(z=0., x=0., lib=np):
      '''
      Returns function f at a given (set of) point(s).
      '''
      f2 = z*0. # allocate and zero
      m = (x**2 + z**2)/2
      f2 = ellipk(m) + ellipe(m)

      return f2

z = df.interpolate(df.Expression('x[0]',degree=cg_degree), V)
x = df.interpolate(df.Expression('x[1]',degree=cg_degree), V)

f1 = df.project(f1(z,x,lib=ufl), V)

plt.colorbar(df.plot(f1))

f2 = df.project(f2(z,x,lib=ufl), V)

plt.colorbar(df.plot(f2))

I define two python functions f1 and f2. In f1 the syntax is the same both for python and ufl and it allows me to project the function on to my FunctionSpace and returns the expected result. However, as f2 is defined in terms of the scipy functions it returns an error when I try to project:

TypeError: ufunc 'ellipk' not supported for the input types

The goal would be to include this function in the weak form of my problem.

4

You should create a UserExpression to do this computation:

5

Hi @dokken thanks for the indication it is indeed what I was looking for. However, I’m encountering a problem with it, when I run the code below, the values for the ellipe_ufl are correct but the nodal values for ellipk_ufl are a string nan. This is not the case if I use ellipk_ufl = ellipk_ufl(degree = 0).

import dolfin as df
from scipy.special import ellipk, ellipe

# Define the mesh and function spaces

mesh_line = df.UnitIntervalMesh(50)
V_line = df.FunctionSpace(mesh_line, 'CG', 1)

k2 = df.interpolate(df.Expression("x[0]", degree = 1), V_line)

class ellipk_ufl(df.UserExpression):
      def set_k2(self, k2):
            self.k2 = k2
      def eval(self, value, x):
            value[0] =  ellipk(self.k2(x))

class ellipe_ufl(df.UserExpression):
      def set_k2(self, k2):
            self.k2 = k2
      def eval(self, value, x):
            value[0] =  ellipe(self.k2(x))

ellipk_ufl = ellipk_ufl(degree = 2)

ellipk_ufl.set_k2(k2)

ellipk_ufl = df.project(ellipk_ufl, V_line)

ellipe_ufl = ellipe_ufl(degree = 2)

ellipe_ufl.set_k2(k2)

ellipe_ufl = df.project(ellipe_ufl, V_line)

print(ellipk_ufl.vector()[:])
print(ellipe_ufl.vector()[:])
6

For the same of clarity, please use distinct variable names in:

as this is quite hard to read as is.

Please note that whenever you use Expression with a given degree, this is interpolated into the appropriate function space. As your function ellipk is note defined at x=1, this causes your interpolation to be nan, making the projection fail. This can be illustrated with:

ellipk(V_line.tabulate_dof_coordinates())

yielding

array([[       inf],
       [3.35414145],
       [3.01611249],
       [2.8207525 ],
       [2.68355141],
       [2.57809211],
       [2.49263532],
       [2.42093296],
       [2.35926355],
       [2.30523174],
       [2.25720533],
       [2.2140225 ],
       [2.17482709],
       [2.13897018],
       [2.10594832],
       [2.07536314],
       [2.04689408],
       [2.02027943],
       [1.99530278],
       [1.97178316],
       [1.94956775],
       [1.92852632],
       [1.90854702],
       [1.88953308],
       [1.87140024],
       [1.85407468],
       [1.83749136],
       [1.82159273],
       [1.80632756],
       [1.79165012],
       [1.77751937],
       [1.76389839],
       [1.7507538 ],
       [1.73805537],
       [1.72577561],
       [1.71388945],
       [1.70237398],
       [1.6912082 ],
       [1.68037282],
       [1.66985009],
       [1.6596236 ],
       [1.64967821],
       [1.63999987],
       [1.63057555],
       [1.62139314],
       [1.61244135],
       [1.60370965],
       [1.59518822],
       [1.58686785],
       [1.57873991],
       [1.57079633]])