This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Assemble a form for integral kernel

+4 votes

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.

asked Feb 25, 2015 by luftiq FEniCS Novice (560 points)
edited Feb 25, 2015 by luftiq

1 Answer

+3 votes
 
Best answer

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.

answered Mar 1, 2015 by MiroK FEniCS Expert (80,920 points)
selected Mar 2, 2015 by luftiq

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.

...