Inner product of grid functions / Access to trial functions

Main Forum Old-Forum - do not post here
1

Hi everybody,

I want to compute the inner product of two grid functions u_i and u_j over a finite element space, but the inner product

InnerProduct(ui,uj)

returns yet another coefficient function instead of the expected scalar
[tex]
\int_{\Omega}u_i(x)u_j(x) \quad .
[/tex]
I was also wondering if I could use a precomputed matrix, containing the inner products of the trial functions, which I tried to construct by

u = fes.TrialFunction() 
v = fes.TrialFunction()
  
M = BilinearForm(fes, symmetric=True)
M += SymbolicBFI(u*v)        
M.Assemble()

Am I missing the point here? Because I thought that’s exactly what a symbolic integrator would do? Is there a way to access the Trial functions by index?

kind regards and thanks in advance
Lukas

2

Hi Lukas,

You can simply call val = Integrate(gfu1*gfu2, mesh)
to compute the inner product.

Best,
Guosheng

3

Hi Guosheng,

thanks for your reply. I also considered this as a possibility, but I am doing this for a large amount of functions u_i and u_j and the Integration is rather slow.
Although this technically solves the issue, I wonder if there is a faster way to do it.
In my understanding, it should be possible to express every Integration by means of an Integration of the trial functions, which has to be computed only once.

4

As you said, just use a BilinearForm:

u, v = V.TnT()
mass = BilinearForm(V)
mass += u*v*dx
mass.Assemble()

tv = ui.vec.CreateVector()

tv.data = mass.mat * ui.vec
ip = InnerProduct(tv, uj.vec)

This way you only need to do one sparse matrix vector multiplication and one InnerProduct of vectors everytime. You just need to assemble the BilinearForm once.

One additional remark:

BilinearForm(..., symmetric=True)

will give you a matrix that is stored symmetrically. This means you use less memory, but the multiplication is a bit slower. If you have to do this for many functions, not using symmetric storage might pay off.

Best,
Lukas

5

Hi Lukas,

yes this does the trick! Thank you very much :slight_smile:
And thanks for the remark, regarding the large amount of functions I use, this may definitely speed up my program.
Accessing the Trial functions by index of the nodes would still be a topic of interest for me, but I can move on with my code now.

Best
Lukas

6

You can access them using u.vec[i] but note that except for p1 H1 these values are not nodal ones, but coefficients for the (high order) basis functions.
You can get these values as a numpy array as well with u.vec.FV().NumPy()