# How to include a known array into a PDE using Expression() Syntax?

I am looking to solve a PDE over a square mesh. In the PDE, there is a function which I’ve called P() and it outputs a Gaussian centred in the middle of the square. But for some cases, I have an array of the values for P() and I don’t know the symbolic formula. How can I replace the function P() with the matrix P_mat in the PDE? A MWE is below.

I am running FEniCS version 2019.1.0 on Docker.

Thank you!

from dolfin import *
import numpy as np

mesh = RectangleMesh(Point(-1,-1), Point(1,1), 50, 50)

V = FunctionSpace(mesh, 'CG', 2)

u = Function(V)
v = TestFunction(V)

P_mat = np.random.rand(50,50) # this is my known array that I want to put into the PDE

def P():
x = Expression('exp((-pow(x[0],2)-pow(x[1],2))/(2*0.01))', degree = 2)
return x

vtkfile = File('MWE/Solution.pvd')

solve(F == 0, u)
vtkfile << u



First you need to explain how the matrix P_mat is supposed to related to the coordinates (x[0], x[1]).
Does the ith row correspond to y = -1+2/50\cdot i, and the jth column to x=-1+2/50\cdot j?

Yes, that is correct. The matrix P_mat is on the coordinates of the mesh, which in this case would be the formulas you pointed out.

Consider the following:

from dolfin import *
import numpy as np

N = 25
mesh = RectangleMesh(Point(-1,-1), Point(1,1), N, N)

V = FunctionSpace(mesh, 'CG', 2)

u = Function(V)
v = TestFunction(V)

P_mat = np.ones((N+1, N+1))
P_mat[:, 0] = 5
P_mat[2, :] += 3

class P(UserExpression):
def __init__(self, data, N,  **kwargs):
super().__init__(kwargs)
self.data = data
self.N = N

def eval(self, value, x):
i = int(self.N/2*(x[1]+1))
j = int(self.N/2*(x[0]+1))
value[0] = self.data[i, j]

def value_shape(self):
return ()

File("P.pvd") << project(P(P_mat, N), V)