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

I/O
1

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

F = dot(grad(v),grad(u))*dx + v*u*dx - v*P()*dx

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

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

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?

3

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.

4

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)

F = dot(grad(v),grad(u))*dx + v*u*dx - v*P(P_mat, N)*dx

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

solve(F == 0, u)
vtkfile << u
3 Likes
5

Yes! Can confirm that this works. Thank you!