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.
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,2)-pow(x,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