Using a spline for nonlinear PDE's coefficient

variational formulation
1

I want to use a spline function in a PDE. An example is the nonlinear Poisson problem given on FEniCS tutorial Solving PDEs in Python - <br> The FEniCS Tutorial Volume I
where instead of the (1+u**2) in q(u), I would like to use a spline (such as a scipy spline interpolation). For example, for a comparison, I was trying

from scipy.interpolate import UnivariateSpline
import numpy as np
un = np.linspace(0,10,20)
q = UnivariateSpline(un,1+un**2,k=2)

def q_analytical(u):
    "Return nonlinear coefficient"
    return 1 + u**2

#The rest of the code is the same as the tutorial example
from fenics import *
# Use SymPy to compute f from the manufactured solution u
import sympy as sym
x, y = sym.symbols('x[0], x[1]')
u = 1 + x + 2*y
f = - sym.diff(q_analytical(u)*sym.diff(u, x), x) - sym.diff(q_analytical(u)*sym.diff(u, y), y)
f = sym.simplify(f)
u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)
print('u =', u_code)
print('f =', f_code)

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression(u_code, degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = Function(V)  # Note: not TrialFunction!
v = TestFunction(V)
f = Expression(f_code, degree=2)
F = q(u)*dot(grad(u), grad(v))*dx - f*v*dx

# Compute solution
solve(F == 0, u, bc)

However, the above code gives an error because u is a ufl object, while the UnivariateSpline expects a float or numpy array. Any idea on how to solve this?

2

@kamensky is an expert in this area. But I’ll try not to embarrass myself with the following response:

Splines are awkward to define symbolically since they may be in C^0. Consider the case where you have repeated knots in the knot vector. The scipy class you’re using is a numerical function, and obviously shouldn’t be used with symbolic quantities as you’ve pointed out. Just browsing the sympy documentation there’s some support for splines. It might not take too much work to develop your own implementation in ufl. Or use sympy's lamdify to generate ufl expressions.

2 Likes