Using FEniCS to generate sampling matrix

For solving inverse problems, we usually need to sample a function at some discrete points. Considering the following problem
$d = \mathcal{F}u + \epsilon,$
with $\mathcal{F} := \mathcal{S} \cdot F$ . Here $F$ represents a partial differential equation or other forward operator, $\mathcal{S}$ is a sampling operator can be defined as
$\mathcal{S}(f) = (f(x_1), f(x_2), \cdots, f(x_{N_{d}}))^{T}.$

Usually, for linear operator $F$ , it can be discretized as a matrix also denoted as $F$ . In the following, we provide a strategy to generate the sampling matrix $\mathcal{S}$ . For finite element method, a function can be expanded by a series of basis functions
$u = \sum_{i=1}^{N}u_{i}e_{i}.$
Then we have
$\left[\begin{array}{c} u(x_1) \\ \vdots \\ u(x_{N_d}) \end{array}\right] = \left[\begin{array}{ccc} e_{1}(x_1) & \cdots & e_{N}(x_1) \\ \vdots & & \vdots \\ e_{1}(x_{N_d}) & \cdots & e_{N}(x_{N_{d}}) \end{array}\right] \left[\begin{array}{c} u_1 \\ \vdots \\ u_N \end{array}\right]$
Hence, the sampling matrix can be obtained by the basis functions evaluated at sampling points. Notice that
$u(x_1) = \sum_{i = 1}^{N}u_{i}e_{i}(x_{1}),$
so if we take $u_1 = 1$ and $u_i = 0$ for $i=2,3,\cdots, N$ , we find
$u(x_1) = e_{1}(x_1).$
Based on this analysis, we can generate a function as follow:

generate a function u in FEniCS;
take each component of u.vector() to be 1 and other components to be 0, and evaluate the function u at the required point.

In the following, we provide a program. The parameter V can be taken as

V = fe.FunctionSpace(mesh, 'P', 1)

The parameter points can be taken as

points = [(0.5,0.5), (0.5, 0.6)]

The program:

import numpy as np
import fenics as fe

def measureMatrix(V, points):
    len_points = len(points)
    u = fe.Function(V)
    u.vector()[:] = 0
    u_vec = u.vector()[:]
    len_u = len(u_vec)
    sampleMatrix = np.zeros((len_points, len_u))
    
    for i in range(len_points):
        for j in range(len_u):
            u_vec[j] = 1
            u.vector()[:] = u_vec
            sampleMatrix[(i, j)] = u(points[i])
            u_vec[:] = 0
    
    return sampleMatrix

The following program illustrates the effectiveness of the above program.

import fenics as fe
import numpy as np

mesh = fe.UnitSquareMesh(10, 10)
V = fe.FunctionSpace(mesh, 'P', 1)

u_D = fe.Expression('1+x[0]*x[0]+2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = fe.DirichletBC(V, u_D, boundary)

u = fe.TrialFunction(V)
v = fe.TestFunction(V)
f = fe.Constant(-6.0)
a = fe.dot(fe.grad(u), fe.grad(v))*fe.dx
L = f*v*fe.dx

u = fe.Function(V)
fe.solve(a == L, u, bc)

points = []
for i in range(100):
    points.append((0.5, i/100.0))
    
S = measureMatrix(V, points)
val1 = S@u.vector()[:]
val2 = [u(i) for i in points]
print(np.linalg.norm(val1-val2))

Using FEniCS to generate sampling matrix

推荐阅读更多精彩内容