I would like to generate noise in a certain way. Assume I have my FEM-space V_h.
Then my noise is to be generated via
\sum_{k=1}^{N_h} \xi_k e_k
where e_k is a basis for the FEM space and \xi_i are normally distributed variables. The issue is to obtain a basis for the FEM space. This method has to be generalizable, so I wonder if there is some “generate basis”-funtion for the FEM-space? There seems to be something like it, https://fenicsproject.org/docs/dolfin/1.6.0/python/programmers-reference/cpp/fem/BasisFunction.html
but I don’t understand the documenation, and more importantly, how do I get this to work with fenics? I need some sort of “integrable” expression, right? Furthermore, I want the covariance matrix R of the noise to be given by
from fenics import *
import numpy as np
# mesh, FunctionSpace
mesh = UnitIntervalMesh(5)
Vh = FunctionSpace(mesh, 'P', 1)
# allocate function u in Vh
u = Function(Vh)
Nh = u.vector()[:].size
xi = np.random.rand(Nh, ) # or any other random array function
# set the entries on the nodes of u to (the random array) xi
u.vector()[:] = xi
# generating mass matrix
z = TestFunction(Vh)
y = TrialFunction(Vh)
m = y*z*dx
Mass = assemble(m)
# Mass.array() gives you the desired matrix
M = Mass.array()
L = np.linalg.cholesky(M)
R = np.transpose(L)
I do not know if the last two lines gives you the covariance matrix R. I supposed that R should be upper triangular? The Cholesky-decomposition gives a lower triangular matrix L such that L L^T = M. Define R := L^T, which is upper triangular, then we have R^T R = L^{T^T} L^T = L L^T = M.
I am sorry, I am new to fenics (not helped by the complete lack of documentation), but I don’t see how that random generation would generate according to my formula above… it is quite important that it takes a sum of normal random coefficients times an eigenbasis function, that is to say
\sum_{k=1}^{N_h} \xi_k e_k
where then e_k is the eigenbasis.
Furhtermore, I want to use this as a RHS in an elliptic equation. I can’t get this to a form that fenics accepts or understand.
Ok. So let u = \sum_{k = 0}^{N_h} u_k e_k, where e_k is your basis of V_h. With u = Function(V) I get exactly such function in FEniCS. Here u.vector()[0] corresponds to u_0, u.vector()[5] to u_5 and so on. With u.vector()[:] I get an array U = [u_0, u_1, \ldots, u_{N_h}]. What I want to do now is to set u_k = \xi_k for k = 0, \ldots, N_h. Thats why I am generating a random-array xi. Now I have to assign the nodal values of u with the of \xi. Translating this is u.vector()[k] = xi[k] for all k. Or, in short notation u.vector()[:] = xi. Here, I assign an array with an array.
Hope that helps. Otherwise I have not understood your question.