Assemble a form for integral kernel

Question

How would one assemble a form so that the result is a matrix

$$K_{ij}=\int_D \int_D C(\boldsymbol{x},\boldsymbol{y}) \phi_i(\boldsymbol{y}) \phi_j(\boldsymbol{x})\, d\boldsymbol{y} d\boldsymbol{x}$$

where $C(x,y)$ is a function and $\phi_i(x)$ are basis for the finite element space.

This is basically related to the solution of the Fredholm integral eigenvalue problem.

Answer 1

Hi, this is just an idea but suppose you set $F_i(x)=\int_{D}C(x, y)\phi_i(y)dy$ so that
$K_{ij}=\int_{D}F_i(x)\phi_j(x)dx$. Now if you replaced each $F_i(x)$ by its interpolant on $V$ then $K_{ij}=F_{ik}M_{kj}$, where $M$ is the mass matrix and $F$ is a matrix of the interpolants. For example for Lagrange elements, computing $F$ would mount to computing
$\int_{D}C(x=x_k, y)\phi_j(y)dy$ where $x_k$ is the coordinates of $k-$th dof. This can be done by assembling a linear form.

Comments

Thanks, this helps. As I am not a FEniCS expert, I am having some trouble implementing this. Here is what I have

import numpy as np
from dolfin import *

mesh = Mesh("boxmesh.xml")
V = FunctionSpace(mesh,"Lagrange",1)
x = SpatialCoordinate(mesh)
gdim = mesh.geometry().dim()
Vdim = V.dim()
dofmap = V.dofmap()
dofs = dofmap.dofs()
dofs_x = dofmap.tabulate_all_coordinates(mesh).reshape((-1,gdim))
kernelfoo = Expression('sigma*sigma*exp(-sqrt((x[0]-y0)*(x[0]-y0) + (x[1]-y1)*(x[1]-y1) + (x[2]-y2)*(x[2]-y2))/l)', sigma = Constant(1.0),l = Constant(1.0), y0 = 1.0, y1 = 1.0, y2 = 1.0)

u = TrialFunction(V)
v = TestFunction(V)

m = inner(u,v)*dx
A = PETScMatrix()
M = PETScMatrix()
assemble(m, tensor=M)

for dof, dof_x in zip (dofs, dofs_x):
    kernelfoo.y0 = dofs_x[dof][0]
    kernelfoo.y1 = dofs_x[dof][1]
    kernelfoo.y2 = dofs_x[dof][2]
    row = np.array([dof], dtype = np.intc)
    fi = inner(kernelfoo,v)*dx
    Fi = assemble(f)

So, basically I am assembling a linear form F_i for each dof. I would then like to assemble these into a matrix and get F. However, I don't know how to manipulate the PETSc matrices...

---

Hi, I did it with uBLAS matrices. See here.